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ABSTRACT 

A  two-scale  analysis  of  the  forced  Rayleigh-Plesset  equation  of  cavitation  bubble 
dynamics  is  performed.  The  problem  of  cavitation  inception  as  it  relates  to  bubble  dynam¬ 
ics  involves  vaporous  cavitation  nucleus  growth  as  it  is  influenced  by  the  pressure 
distribution  on  a  submerged  body.  This  brings  into  prominence  two  widely  varying  time 
scales.  The  "laboratory  time"  is  characterized  by  the  bubble’s  travel  over  the  body,  while 
the  "bubble  time"  is  characterized  by  the  very  high  natural  frequency  oscillations  of  the 
individual  bubble.  The  laboratory  time  is  expected  to  be  much  longer  than  the  bubble  time; 
thus  they  can  be  related  by  a  very  small  parameter,  e. 

Using  these  two  time  scales,  a  perturbation  expansion  is  performed  on  the  forced 
Rayleigh-Plesset  equation  and  its  initial  conditions  up  to  the  third  order  in  e.  The  resulting 
zero  and  first-order  equations  are  solved,  subject  to  these  solutions  being  independent  of  the 
laboratory  time.  In  this  case  the  integrability  condition  for  each  step  is  thereby  identically 
satisfied. 

The  zero-order  equation  is  an  autonomous  nonlinear  second-order  differential  equa¬ 
tion.  An  approximate  closed-form,  small  oscillation  solution  is  found  and  a  numerical 
solution  is  found  for  oscillations  of  all  amplitudes.  Both  forms  of  the  solution  are  found  to 
be  periodic  and  a  function  of  bubble  time  only. 

The  first-order  equation  is  a  nonhomogeneous  linear  second-order  differential 
equation  with  periodic  coefficients.  Its  homogeneous  form  is  of  the  class  of  Hill's  equa¬ 
tions,  and  can  be  treated  using  Floquet  theory.  The  complementary  solution  in  normal  or 
Floquet  form  is  found  numerically.  Variation  of  parameters  is  used  to  find  the  particular 
integral.  Again,  the  constants  fa*  the  first-order  solution  are  independent  of  the  laboratory 
time  and  only  the  bubble  time  appears  explicitly  in  the  complementary  function. The  labora¬ 
tory  time  enters  the  solution  via  the  particular  integral  because  the  forcing  function  depends 
cm  the  laboratory  time. 
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The  dynamic  stability  of  these  solutions  was  investigated.  It  was  found  that  the 
present  solutions  are  unstable  and  they  are  not  suitable  for  inception  calculations.  Further 
study  of  the  method  of  solution  determined  that  a  more  general  solution  than  the  present  one 
might  be  found  which  will  still  satisfy  the  integrability  conditions  at  each  step.  Some  basic 
zero-order  equations  have  been  formulated  which  may  permit  such  a  numerical  solution  to 
be  found.  No  algorithm  for  the  solution  of  these  equations  has  been  developed. 

This  report  is  a  revision  of  the  first  author's  thesis  in  Aerospace  Engineering,  which 
was  submitted  in  partial  fulfillment  of  the  requirements  for  the  degree  of  Master  of  Science, 
December  1990. 
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CHAPTER  1 
INTRODUCTION 


1.1  Background 

In  a  flow  about  a  submerged  body,  the  fluid  is  known  to  contain  very  small  nuclei 
containing  air  and  water  vapor.  These  nuclei  translate  downstream  at  some  velocity  close 
to  that  of  the  free  stream.  Some  of  these  will  come  in  contact  with  the  body  and  enter  its 
boundary  layer.  Submerged  bodies,  such  as  circular  cylinders,  hydrofoils  or  hemispherical 
headforms,  have  pressure  distributions  that  vary  along  their  arc  lengths.  At  some  point  on 
the  body  the  local  static  pressure  of  the  liquid  will  equal  that  of  the  liquid  vapor  pressure 
after  which  there  will  be  a  region  in  which  the  static  pressure  is  less  than  vapor  pressure.  It 
is  in  this  region  that  the  nuclei  can  experience  vaporous  bubble  growth.  Therefore  the 
critical  point  at  which  the  static  pressure  equals  the  vapor  pressure  separates  two  possible 
regions  of  bubble  behavior  along  the  body. 

If  the  nucleus  enters  the  boundary  layer  of  the  body  upstream  of  this  critical  point,  it 
will  travel  along  the  body  and  act  as  a  flaccid  bubble.  In  this  upstream  region  the  bubble 
will  grow,  mainly  due  to  internal  gas  pressure  changes  caused  by  the  changing  liquid  static 
pressure,  until  it  sees  a  local  static  equal  to  vapor  pressure.  After  this  point  on  the  body  the 
local  static  pressure  is  less  than  the  vapor  pressure  within  the  bubble  and  the  conditions  are 
favorable  for  vaporous  bubble  growth.  If  the  nucleus  either  enters  the  boundary  layer  of 
the  body  in  this  region  or  is  conveyed  from  the  flaccid  bubble  region  into  the  favorable 
region,  it  will  grow  vaporously  to  a  maximum  radius  and  then  collapse  as  the  local  static 
pressure  rises. 

If  there  is  no  separation  on  the  body  the  bubble  will  continue  to  collapse  rapidly  and 
violently.  For  sufficiently  low  Reynolds  numbers,  however,  some  bodies  will  have  lami¬ 
nar  separation  regions.  If  the  separation  bubble  is  short,  it  seems  possible  that  collapse 
may  not  occur,  and  the  bubble  will  come  to  rest  within  the  laminar  separation  bubble  where 
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it  can  grow  by  gas  diffusion  [1],  [13].  This  diffusive  growth  continues  until  the  bubble 
diameter  reaches  the  height  of  the  separation  zone.  Then  the  bubble  interacts  with  the  free 
shear  layer  and  is  translated  downstream.  If  the  flow  reattaches,  the  turbulent  shear  in  this 
region  breaks  the  bubble  into  froth,  creating  a  narrow  ring  of  visible  cavitation.  The  point 
at  which  this  cavitation  becomes  visible  is  called  cavitation  inception.  This  type  is  known 
as  bubble-ring  cavitation,  which  is  a  form  of  attached  cavitation  controlled  primarily  by 
laminar  boundary  layer  separation. 

The  governing  equation  which  describes  isothermal  cavitation  vapor  bubble  growth 
or  collapse  is  the  Rayleigh-Plesset  equation  [16].  It  describes  the  growth  of  a  spherical 
bubble  in  the  region  of  a  body  where  the  local  static  pressure  is  less  than  the  liquid  vapor 
pressure.  It  also  applies  at  least  in  the  initial  stages  of  collapse  when  the  static  pressure  rises 
down  stream  of  the  minimum  pressure  point  on  the  body.  An  examination  of  this  equation 
may  provide  direct  insight  into  the  hydrodynamically  forced  growth  of  the  nuclei  leading 
to  cavitation  inception.  From  this  knowledge  it  is  hoped  that  a  relationship  between  the 
parameters  of  the  flow  and  submerged  body  with  respect  to  critical  flow  conditions  for 
cavitation  onset  might  be  found. 

1.2  Previous  Investigations 

The  governing  Rayleigh-Plesset  equation1  describing  isothermal  vaporous  growth 
and  collapse  of  a  spherical  bubble  of  radius,  R,  as  a  function  of  time,  t,  in  a  perfect  fluid,  is 
written  as 

r,  d2R  ^  3  fdRf  x  _1_  .  tj/.n 
R  dt2  2  [dt  J  "  R3  '  R  +  F(t)  •  (U) 

Equation  (1.1)  assumes  that  the  static  pressure  is  responsible  fen*  driving  the  bubble 
growth,  and  is  contained  in  the  time  dependent  forcing  pressure ,  F(t). 


1See  Appendix  D  for  more  dedails. 
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that  the  forcing  function  pulse  acting  on  a  nucleus  is  characterized  by  two  widely  varying 
time  scales.  One  dimensionless  scale,  called  "laboratory  dme,"  t,,  is  governed  by  the  time 

it  takes  for  the  growing  nucleus  to  move  through  the  low  pressure  region  in  the  boundary 
layer  which  promotes  vaporous  growth.  The  other  dimensionless  time  scale,  called 
"bubble  time,"  x,  is  associated  with  the  very  high  natural  frequency  of  the  very  small 
nuclei  and  depends  on  the  period  of  natural  vibration  of  the  a  typical  microbubble.  This 
guarantees  a  very  rapid  response  of  the  bubble  to  external  pressure  changes.  For  every 
period  of  "slow"  laboratory  time  there  are  many  periods  of  the  "fast"  bubble  time.  There¬ 
fore,  the  two  time  scales  can  be  related  by  a  very  small  parameter,  e,  which  is  simply  the 
ratio  of  ts  to  x. 

Using  this  idea.  Baker  expanded  the  forced  Rayleigh-Plesset  equation  and  its 
initial  conditions  in  orders  of  the  small  parameter  e.  He  then  found  an  approximate 
solution  to  the  zero-order  equation  in  terms  of  elliptic  integrals  and  functions.  Although 
this  analysis  was  quite  complex  and  involved  some  approximations,  it  compared  well 
with  a  strict  numerical  solution  of  the  zero-order  equation.  Also,  it  set  up  a  basis  for 
the  analysis  of  the  first-order  equation  of  the  expansion  which  can  be  combined  with 
the  zero-order  findings  to  approximate  better  the  total  solution. 

1.3  Objective  of  this  Investigation 

The  main  objective  of  this  investigation  is  to  find  a  solution  of  the  forced  Rayleigh- 
Plesset  equation  (1.1).  In  general  this  can  be  solved  numerically  using  standard  computer 
algorithms  which  do  not  consider  the  aforementioned  time  scales.  Use  of  this  type  of 
numerical  approach  makes  it  extremely  difficult  to  isolate  any  critical  flow  parameters  or 
conditions  between  them  which  can  lead  to  cavitation  inception. 

litis  investigation  extends  Baker's  method  of  solution  [2].  It  involves  a  two-scale 
perturbation  expansion  based  on  the  small  parameter  E,  which  relates  the  two  time  scales 
present  in  the  problem,  t,  and  x.  Application  of  this  two-scale  procedure  produces  a  series 
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of  linear  or  nonlinear  partial  differential  equations  for  the  various  powers  in  e.  This  type  of 
mathematical  technique  is  well  known  [11]  and  it  is  hoped  that  it  will  work  well  with  the 
very  nonlinear  equation  (1.1). 

The  first  part  of  the  present  investigation  considers  the  forcing  function.  This  comes 
from  the  pressure  distribution  across  a  particular  body.  Baker  used  a  hemispherical  head- 
form  so  that  the  pressure  distribution  may  be  correlated  with  much  experimental  data.  This 
analysis  will  attempt  to  use  a  circular  cylinder  as  the  governing  body,  in  part  to  see  if  the 
analysis  is  valid  for  different  cases.  Supercritical  flow  over  a  circular  cylinder  exhibits  a 
pressure  distribution  showing  a  laminar  bubble  just  upstream  of  the  turbulent  separation 
[3],  [4].  Thus,  experimental  pressure  data  for  a  circular  cylinder  will  be  used  to  derive  a 
forcing  function  for  the  problem. 

Next,  the  two-scale  dynamical  equations  will  be  derived.  Although  Baker's  solu¬ 
tion  of  the  zero-order  equation  gave  some  new  approximations  for  the  time  of  growth  in 
terms  of  the  radius,  it  can  not  be  inverted  explicitly  to  give  the  radius  as  a  function  of  time. 
An  inversion  is  needed  in  order  to  obtain  the  first-order  solution.  Thus,  in  the  present 
analysis,  we  will  rely  on  numerical  methods,  supplementing  them  with  small  oscillation 
analysis  which  leads  to  periodic  solutions  requiring  no  inversion.  Baker  also  assumed  that 
the  zero-order  solution  was  independent  of  the  laboratory  time,  ts.  This  investigation  will 

also  attempt  to  verify  this  fact  or  at  least  to  see  what  further  research  is  needed  if  conditions 
can  be  found  which  show  that  other  zero-order  solutions  may  be  possible. 

Lastly,  the  first-order  equation  is  analyzed  for  the  case  of  periodic  zero-order 
solutions  and  a  numerical  first-order  solution  is  obtained  using  the  Floquet  theory.  First  the 
complementary  function  is  found  and  its  dynamic  and  numerical  stability  is  examined.  Then 
the  Particular  integral  is  found  using  variation  of  parameters.  The  entire  process  is  trans¬ 
lated  into  Fortran  and  a  numerical  example  of  the  complete  solution  is  given. 
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CHAPTER  2 

BACKGROUND  AND  FORMULATION 


2.1  Development  of  the  Time  Scales 

As  mentioned  from  the  outset,  this  problem  involves  two  time  scales.  The  first  of  these  is 
associated  with  the  time  it  takes  for  the  growing  nucleus  to  move  through  the  low  pressure 
region  on  the  submerged  body,  and  the  second  is  associated  with  the  very  high  frequency 
response  time  of  the  nucleus  to  external  pressure  changes. 

It  is  desirable  to  have  the  two  time  scales  in  dimensionless  form.  By  first  defining  a 
characteristic  laboratory  time,  D/Vq,  one  can  define  the  dimensionless  laboratory  time  as 


tVo 

D 


(2.1) 


where  t  is  real  time,  D  is  the  characteristic  diameter  or  other  length  of  the  body  and  V0  is  the 
free  stream  velocity.  Now  we  define  a  characteristic  time,  -yj  p  R^/2a  ,  involving  only 


the  physical  properties  of  the  bubble  and  liquid.  This  characteristic  time  enables  one  to 
define  a  dimensionless  bubble  time  as 


x 


(2.2) 


where  p  is  the  liquid  density  and  o  is  the  liquid  surface  tension.  Introduction  of  a  small 
parameter  e  as  the  ratio  of  the  characteristic  bubble  time  to  the  characteristic  laboratory  time 


yields: 


Vo  A  (2-3) 

e  =  D/V0  =  D  V  2o  ' 

It  should  also  be  noted  that  such  a  definition  also  gives  the  ratio  £  in  terms  of  the  nondi- 
mensional  terms  defined  above,  i.e.. 


e  =  tj/t. 


(2.4) 
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A  relative  comparison  of  the  "slow"  laboratory  time  scale  and  the  "fast"  bubble  time 
scale  shows  that  the  ratio,  x/tj ,  can  have  a  magnitude  of  about  103  or  less.  This  says  that  x 
is  of  very  short  duration  compared  to  the  time  scale  tj,  hence  there  are  many  periods  of  x 
for  one  period  of  tj.  It  follows  from  equation  (2.3)  that  e  should  be  a  small  number,  often 
of  the  order  of  10'3  .  These  are  the  same  two  time  scales  used  by  Baker  [2]  in  his  work. 

As  he  states,  the  fast  bubble  time  characterizes  the  high  frequency  oscillations  of  the 
individual  bubbles,  whereas  the  slow  laboratory  time  characterizes  the  forcing  function 
across  the  governing  body. 


2.2  The  Forcing  Function 

In  a  discussion  of  the  translation  of  a  bubble  over  a  submerged  body,  what  must  be 
brought  to  light  is  the  driving  force  which  creates  the  environment  for  subsequent  vaporous 
bubble  growth  and  the  start  of  collapse.  Discussion  has  been  made  of  the  fact  that  for  va¬ 
porous  growth  to  occur  there  must  be  a  region  on  the  body  where  the  static  pressure  of  the 
flow  is  less  than  the  vapor  pressure  of  the  liquid.  This  region  starts  on  the  body  at  a  critical 
point  where  the  static  pressure  equals  that  of  the  vapor  pressure  of  the  liquid.  Downstream 
of  this  critical  point  in  the  growth  region,  the  pressure  is  characterized  by 


Pv  >  P  *  Pmin  .  (2.5) 

where  pv  is  the  vapor  pressure  of  the  liquid  and  p,^  is  the  minimum  static  pressure  on  the 

body. 

Defining  a  pressure  coefficient  as  a  function  of  arc  length  along  the  body  and  the 
cavitation  number  respectively  as 


Cp(s)  = 


p(s) ■  Po 


(2.6) 


and 
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K  = 


PO  '  Pv 


(2.7) 


one  can  write  the  relationship  in  equation  (2.5)  as 

K  >  *  *-p(s)  ^  *  (Cp)min  •  (2.8) 

Here,  p0  is  the  free  stream  static  pressure.  Equation  (2.8)  suggests  that,  depending  on  the 

value  of  K,  the  negative  of  the  pressure  coefficient  can  be  used  to  measure  the  pressure 
force  on  the  bubble  that  causes  vaporous  growth.  If  Cp  =  -  K,  the  pressure  force  would  be 

zero.  Thus  a  positive  forcing  function  which  acts  on  a  typical  bubble  can  be  defined  as  the 
negative  difference  in  these  two  quantities  or 

F(et)  =  -Cp(s)  -  K  =  -  [Cp(s)  +  K]  .  (2.9) 

It  should  be  noted  in  this  particular  case  that  the  region  of  vaporous  growth  and 
collapse  extends  along  the  arc  of  the  body  only  as  far  as  the  separation  point  of  the 
laminar  region.  Once  the  bubble  enters  the  separation  region,  it  is  assumed  that  the 
process  of  vaporous  growth  ends  and  any  growth  that  continues  is  due  to  gas  diffusion 
from  the  liquid  into  the  bubble.  Figure  2.1  shows  a  schematic  of  the  pressure 
distribution  across  a  hemispherical  headform  borrowed  from  Baker  [2].  It  gives  an 
example  of  the  different  states  a  bubble  might  see  if  it  traveled  across  the  entire  body 
including  the  flaccid  region,  the  dividing  point  where  Cp(s)  =  -  K,  the  vaporous 

growth  region  and  the  laminar  separation  point 

The  forcing  function  is  written  as  F  =  F(t^  =  F(et)  due  to  the  fact  that  it  acts 

over  the  low  pressure  or  vaporous  growth  region  of  the  body  which  is  characterized  by 
the  laboratory  time  scale.  As  an  aside,  intuition  says  that  the  forcing  function  depends 
on  a  parameter  that  contains  e,  and  it  is  shown  later  that  this  is  the  case.  Since  the 
pressure  coefficient  is  given  as  a  function  of  arc  length  on  the  body,  one  must  convert 
from  the  arc  length  to  the  laboratory  time. 


AVERAGED  PRESSURE  COEFFICIENT, 
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Figure  2. 1  Schematic  Diagram  of  the  States  as  Seen  by  a  Bubble  Traveling 
Along  a  Hemispherical  Headfonn,  after  Baker  [2]. 
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First,  a  dimensionless  arc  length  is  defined  as 

2S 

s  =  D  (2.10) 

where  S  is  the  dimensional  arc  length  of  the  body.  Assuming  the  boundary  layer  to  be  a 
vortex  sheet,  one  can  approximate  its  overall  translational  velocity  as  being  one  half  of  the 
local  flow  speed  at  the  edge  of  the  boundary  layer,  or 


v(s)  =  -yWl  -  CD(s)  , 


(2.11) 


which  then  equals  the  convective  speed  of  a  typical  bubble  in  the  boundary  layer.  In  terms 
of  the  body  length  this  convective  speed  is 


,  .  dS  dS  ds 

v(s>  ■  S  ■  5  • 


(2.12) 


Substituting  the  derivative  of  S  from  (2.10)  and  v(s)  from  equation  (2.11)  into  equation 
(2.12)  and  solving  for  ds/dt  one  has 


|  -  tp/l  -  C„(s)  . 


(2.13) 


Recalling  equation  (2.2)  which  relates  the  bubble  time  with  real  time,  one  can  differentiate 
that  equation  to  get 


dx  = 


(2.14) 


Substituting  equation  (2.4)  for  e  and  (2.13)  into  equation  (2.14)  and  integrating  over  s,  one 
finds  that  an  equation  relating  the  dimensionless  arc  length  to  the  laboratory  time  is 


t,  =  ex  = 


=  .vr 


Cp(t> 


-,  i  1, 2, ....  n  . 


(2.15) 


In  this  case  the  suffix  i  corresponds  to  each  pressure  coefficient  vs.  its  arc  length  datum 
point,  Sj.  This  is  the  same  transformation  used  by  Baker  [2]. 
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What  is  now  needed  is  the  pressure  coefficient  data  that  are  appropriate  to  any 
particular  problem.  The  data  must  be  either  given  as  a  function  of  arc  length  along  the  body 
or  be  readily  convertible  to  such.  We  shall  ask  that  the  particular  body  from  which  the  data 
are  taken  contain  the  necessary  laminar  separation  region  although  its  presence  or  absence 
will  not  preclude  vaporous  bubble  growth.  Baker  used  experimental  Cp  values  given  for  a 

hemispherical  headform.  He  determined  that  for  cavitation  numbers  ranging  from  K  =  0.6 
to  K  =  0.7,  the  laminar  separation  points  corresponded  to  positions  of  s  =  0.813  and 
s  =  0.731,  respectively.  He  obtained  very  good  results  using  these  data. 

In  order  to  check  the  range  of  flow  parameters,  such  as  the  cavitation  number,  this 
analysis  will  attempt  to  use  data  from  a  body  containing  a  lower  (Cp)min,  and  thus  be 

compatible  with  higher  values  of  K.  A  circular  cylinder  was  chosen  because  even  at  critical 
Reynolds  numbers  it  exhibits  laminar  separation,  and  there  exist  some  inception  data  [3]  for 
this  body  although  the  inception  is  usually  defined  by  the  first  appearance  of  cavitation 
bubbles  within  vortex  cores  in  the  separated  shear  layer .  Reference  [3]  reports  observa¬ 
tions  of  cavitation  in  the  laminar  bubble.  The  pressure  distribution  used  in  this  analysis 
comes  from  an  NACA  Technical  Note  by  Gowen  and  Perkins  [4].  Their  graphical  data 
showed  a  Cplmin  of  about  -  2.6  for  critical  Reynolds  numbers  and  at  very  low  Mach 

numbers.  The  tabulated  values  used  are  given  in  appendix  B,  table  B.l,  and  the  smoothed 
data  as  a  function  of  arc  length  interpolated  from  the  table  are  shown  in  figure  2.2. 

Although  the  original  values  were  simply  read  off  of  a  graph,  the  fact  that  they  exhibit  the 
proper  trend  is  more  important  than  their  accuracy  since  this  analysis  is  somewhat 
qualitative. 

The  location  of  the  laminar  separation  point  is  needed  for  the  circular  cylinder  if  (Hie  is  to 
know  when  the  analysis  is  no  longer  valid.  Applying  Thwaites'  method  and  using  a 
computer  code  given  by  Moran  [5],  we  verified  that  in  water  at  a  free  stream  velocity  of 
V0  *  40  ft/s,  a  circular  cylinder  exhibits  laminar  separation.  Using  this  code,  and 

comparing  its  results  with  the  best  estimates  from  Schlichting  [6]  and  Goldstein  [7],  we 
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Pressure  distribution  on  a  circular  cylinder. 
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Figure  2.2  Experimental  and  Smooth  Pressure  Coefficient  Data  Versus  Dimensionless  Arc 
Length  on  a  Circular  Cylinder  at  Supercritical  Reynolds  Numbers 
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found  separation  to  occur  at  approximately  1 10  degrees  from  the  forward  stagnation  point 
or  at  a  dimensionless  arc  length  of  s  =  1.919.  Knowing  this,  one  finds  a  range  of  cavitation 
numbers  that  will  allow  favorable  bubble  growth  and  collapse  to  be  from  K  =  2. 1  to 
K  =  2.4. 

A  Fortran  code  to  determine  the  forcing  function  as  a  function  of  the  laboratory  time 
was  developed  and  is  contained  in  the  complete  code  given  in  appendix  A.  It  begins  by 
curve  fitting  the  given  Cp  vs  s  data  using  a  cubic  interpolation  scheme,  in  order  to  produce 

evenly  spaced  values  so  that  the  integral  in  equation  (2.15)  can  be  evaluated  using  Simp¬ 
son’s  rule.  This  calculation  correlates  the  arc  length  (radians)  with  the  dimensionless 
laboratory  time.  Care  must  be  used  when  evaluating  the  integral  at  the  stagnation  point 
because  of  the  square  root  singularity  involved.  A  parabolic  approximation  is  used  over  the 
first  three  points,  including  the  singularity,  to  take  care  of  this  problem.  The  corresponding 
value  of  the  forcing  function  is  found  by  applying  the  chosen  cavitation  number  to  equation 
(2.9)  for  the  interpolated  Cp  values.  This  step  yields  the  forcing  function  as  a  function  of 

laboratory  time. 

To  give  us  a  feel  for  a  likely  laboratory-time  durations ,  figure  2.3  shows  the 
pressure  coefficient  as  a  function  of  %  as  calculated  from  equation  (2.15).  Figure  2.4  is  a 

schematic  diagram  of  a  typical  forcing  function  showing  its  vaporous  growth  and  collapse 

regions.  The  beginning  of  the  positive  growth  region  is  located  at  the  intersection  of  the 
horizontal  line  at  Cp  =  -  K  and  and  it  is  taken  to  be  the  origin  of  forcing  function  coordi¬ 
nates.  Since  this  is  the  critical  point  where  vaporous  growth  begins  and  the  start  of  the 
region  over  which  this  investigation  takes  place,  this  is  the  point  defined  to  be  where  the 
two  time  scales  are  equal  to  zero.  It  is  convenient  then  to  shift  the  coordinate  axes  to  this 
point  and  have  it  be  die  origin.  Depending  on  the  value  of  K  chosen,  this  shift  will  occur 
along  a  vector  as  shown  in  figure  2.4.  The  computer  code  developed  also  finds  this  new 
origin  and  so  it  automatically  shifts  the  time  axis  so  that  it  begins  at  zero.  Figure  2.5  below 
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Pressure  distribution  on  circular  cylinder. 


Figure  2.3  Pressure  Distribution  on  a  Circular  Cylinder  as  a 
Function  of  the  Labortatory  111116. 


FORCING  FUNCTION. 
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Figure  2.4  Schematic  Diagram  of  a  Typical  Forcing  Function  Showing 
Regions  of  Vaporous  Growth  and  Collapse. 


Forcing  function,  F(ts)  =  -  (Cp+  K). 
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shows  the  shifted  forcing  function  for  a  circular  cylinder  at  values  of  K  =  2.1  and  K  =  2.4 
as  a  function  of  the  actual  time.  The  actual  time  can  be  divided  by  the  ratio  of  D/V0  in  order 

to  get  F  in  terms  of  the  laboratory  time,  which  is  then  used  as  an  input  to  the  bubble- 
dynamic  analysis. 

2.3  Flaccid  Bubble  Analysis 

As  described  by  Baker  [2],  the  bubble  grows  as  a  flaccid  balloon  before  it  reaches 
the  point  where  Cp  =  -K.  It  experiences  successive  states  of  isothermal  equilibrium  as  its 
size  changes  from  the  time  it  reaches  a  point  near  the  stagnation  point  on  a  submerged 
body,  and  it  travels  within  the  boundary  layer  before  it  reaches  the  region  favorable  to 
vaporous  growth.  The  bubble  experiences  no  vaporous  growth  during  this  first  phase  of  its 
motion  on  the  body.  As  sketched  here,  suchbubbles  originate  from  free-stream  micro¬ 
bubbles  of  some  radius  r  at  the  free-stream  static  pressure  p0. 

The  term,  2a/r,  gives  the  pressure  jump  across  the  bubble 
surface.  The  internal  pressure  is  the  sum  of  vapour  and  air 
pressures  as  shown. 

The  flaccid  bubble  kinematics  relates  the  bubble's  radi¬ 
al  motion  to  its  movement  along  the  dimensionless  arc  length, 
s,  of  the  headfoim.  Depending  on  the  size  of  the  submerged 
body,  this  length  is  the  dimensional  arc  length  of  the  head- 

form,  S,  divided  by  some  characteristic  body  length,  D.  The  rate  of  change  of  pressure 
experienced  by  the  bubble  is  therefore, 

dCp  dCp  ds 

A  =  ~5s  dt  *  (2-16) 


A  Free-Stream  Nucleus. 


This  form  is  sought  because  the  experimental  pressure  distribution  is  given  in  terms  of  s. 
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The  kinematics  of  the  previous  section,  as  embodied  in  equation  (2.13),  provides 
an  expression  for  ds/dt.  Substituting  equation  (2.13)  into  equation  (2.16),  one  gets  the  rate 
of  change  of  pressure  experienced  by  a  nucleus  as  it  travels  across  a  body  to  be 


dt  D  *  1  ds  • 


(2.17) 

This  is  all  the  kinematics  that  are  needed.  Next  we  use  Boyle's  law  for  an  isothermal  flaccid 
bubble  and  the  force  (pressure)  balance  across  the  bubble  wall.  Then  we  combine  these 
three  results  to  achieve  expressions  for  the  bubble's  flaccid  growth  characteristics. 

The  bubble  contains  both  air  and  water  vapor.  Therefore  inside  the  bubble  the  total 
pressure  is 


Pi  =  Pa  +  Pv  •  (2.18) 

Only  the  air  obeys  Boyle's  law,  which  requires  that 

P«R3(t)  =  P«oRo  *  (2.19) 

where  p^  =  partial  pressure  of  air  inside  the  free  stream  nucleus  at  R  =  Rq.  A  dimen¬ 
sionless  bubble  radius  is  now  defined  using  the  free  stream  nucleus  radius,  Rq,  such  that 


r(t)  =  R(t)/Ro  .  (2.20) 

It  is  noted  that  p^  is  a  constant,  and  in  the  free  stream  when  r  =  1,  pa  =  p^.  Now  if  p(t) 


is  the  static  pressure  outside  the  bubble  and  a  is  the  surface  tension,  on  the  bubble  wall  the 
pressure  forces  are  in  balance  and 


2a  ,  . . 

Pi  =  +  P(S)  • 

So  for  a  bubble  in  the  free  stream,  as  in  the  illustrative  sketch  above, 

2a 

P*0  +  Pv  =  Rq  +  Po  » 


(2.21) 


(2.22) 


and  for  a  bubble  inside  the  boundary  layer 
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Pa  +  Pv  =  If)  +  p(s)  •  (2.23) 

Solving  (2.22)  and  (2.23)  for  pa  and  respectively  and  substituting  them  into  Boyle's 


law,  equation  (2.19),  one  finds  the  result. 


+  p(s)  -  p 


,]r3(o  =  [i 


+  Po  -  P 


,]rJ  . 


(2.24) 


Bv  substituting  (2.6)  and  (2.7)  for  the  various  pressure  differences  into  (2.24)  and  using 
the  nondimensional  bubble  radius,  one  finds  that 

R^i)  +  5PV0<Cp  +  K)  *,)  =  +  Ip v^K  .  (2.25) 

Multiplying  through  by  Rq/2  o  and  defining  a  radial  Weber  number  as 


(2.26) 


one  obtains  a  cubic  equation  for  the  dimensionless  flaccid  bubble  radius  in  the  form 


iw,[Cp(t)  +  K]  r3(t)  +  m  -  [l  +  ’jKWr]  =  0  .  (2.27) 


This  equation  can  be  solved  exactly  and  its  derivative  taken  to  obtain  expressions  for  the 
flaccid  bubble  radius  and  growth  rate  as  functions  of  real  time,  t.  Because  of  the  fact  that 
the  Rayleigh-Plesset  equation  governs  the  growth  of  the  bubble  after  the  point  of  vaporous 
growth,  this  expression  and  its  derivative  can  be  used  to  fonn  initial  conditions  at  the 
critical  point  where  vaporous  growth  begins. 

In  the  special  case  when  Cp  =  0,  equation  (2.27)  can  be  factored  as  follows: 


[K  W  1 

—L(T2  +  r+  i)  +  r+  1J  (r  - 1)  =  0. 

The  only  root  of  physical  interest  is  r  =  1,  as  should  be  the  case  for  a  nucleus  in  the  free 
stream.  Another  special  case  occurs  when  Cp  =  -K.  Then  the  cubic  term  is  lost,  and  we 
find  that  the  two  remaining  roots  are 
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The  positive  root  is  a  useful  physical  quantity  in  the  dynamical  solution  of  Chapter  3.  We 
note  here  that  it  pertains  to  those  values  of  Cp  and  K  which  mark  the  origin  of  coordinates 
in  figures  2.4  and  2.5.  That  is,  for  any  value  of  K  <  -Cplmi„,  vaporous  growth  first 
becomes  possible;  and  from  equation  (2.20),  it  gives  the  ratio  of  microbubble  size  at  that 
initial  condition  to  a  typical  free-stream  nucleus  size. 


2.3.1  Formulation  of  the  Flaccid  Bubble  Radius.  r(t~) 

The  cubic  flaccid  bubble  equation  can  be  solved  exactly.  Because  of  the  fact  that 
the  non-dimensional  bubble  radius  is  very  close  to  unity,  one  can  easily  find  a  solution 
using  a  simple  perturbation  of  the  form 


r  =  1  -  x  ,  (2.28) 

where  0  <  x  «  1.  Substituting  this  into  (2.27)  and  neglecting  terms  of  order  x3  and 
higher,  one  sees  that  a  quadratic  equation  for  x  remains  and  it  is 

(3 A  +  l)x2  -  (3 A  +  2)x  +  |CpWr  =  0  ,  (2  29) 

where  A  =  ^(Cp  +  K)  Wr .  The  two  roots  of  x  are  found  to  be 


x 


(3A  +  2) 

2  (3  A  +  1)  1 


1  -  CpWr 


(3  A  +  1) 
(3A  +  2)2 


(2.30) 


The  proper  root  must  be  decided  upon  for  this  particular  case.  Consider  the  case  in  which 


Cp  -*  -  K;  then  A  -*  0.  Due  to  the  conditions  imposed  upon  x  by  equation  (2.28),  the 


negative  root  must  be  chosen  to  satisfy  0  <  x  «  1.  The  value  of  x  is  therefore 


lim  x  =  1  -  -v/ 1  +  jKW.. 

A  -  0  V  4  1 

(2.31) 


For  ease  of  writing  a  new  variable  is  introduced,  namely; 
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q-|KWr-  (2.32) 

In  the  limit  as  Cp  =  -  K,  equation  (2.31)  gives  x  =  0  corresponding  to  a  nucleus  of  r  =  1 

as  it  is  in  the  free  stream.  The  same  result  is  found  if  one  applies  this  condition  to  the  cubic 
flaccid  bubble  equation.  Setting  Cp  =  -  K  =  0  in  equation  (2.27),  we  find 

r2  -  1  =  0  ,  (2.33) 

which  has  a  root  at  r  =  1.  For  Cp  =  -  K  *  0  and  neglecting  the  cubic  term  as  before,  the 

equation  for  r  is  written  as 


r2  -  1  -  q  =  0.  (2.34) 

Equation  (2.34)  yields  the  same  result  as  one  gets  by  substituting  (2.31)  into  (2.28), 
nameley: 


r  =  VI  +  q  . 

Substitution  of  (2.30)  into  (2.28),  with  the  expression  for  A  in  mind,  yields 


_  ,  1  (3A  +  2) 

'w  “  1 "  2(3A  +  1) 


1  ±  a/  I  -  C  W 

P  r(3A  +  2)2 


(2.35) 


(2.36) 


This  then  is  the  flaccid  bubble  radius  as  a  function  of  actual  time,  in  terms  of  the  flow 
parameters  Cp,  K,  and  Wr.  It  is  valid  only  in  the  region  from  the  stagnation  point  up  to  the 


initial  point  of  vaporous  growth.  This  is  similar  to  the  form  derived  by  Baker  [2],  and  as 
was  shown  by  him,  shows  very  good  agreement  with  experimental  results  of  nucleus 
measurements  in  this  region. 


2.3.2  Calculation  of  the  Initial  Radius.  rfQl 

For  the  vaporous  growth  phase  the  initial  point  at  which  the  dimensionless  labora¬ 
tory  time  and  the  dimensionless  bubble  time  are  measured  is  the  point  on  the  body  where 
Cp  =  -  K.  Since  vaporous  growth  is  preceded  by  the  flaccid  bubble  region,  the  initial 

conditions  for  vaporous  growth  from  the  flaccid  bubble  equations  should  be  evaluated  at 
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the  point  on  the  body  surface  where  this  Cp  value  is  found.  Recalling  that  equation  (2.35) 
was  derived  at  the  point  where  Cp  =  -  K  and  if  t  =  to  at  this  point,  it  can  be  written  as 

r(to)  =  VT  +  q  .  (2.37) 

This  then  is  the  initial  bubble  radius  at  t  =  to  and  tj  =  t  =  0,  corresponding  to  the  point 
<mi  the  body  where  Cp  =  -  K.  Using  this  result  one  can  define  the  dimensionless  parame¬ 


ter  for  any  real  time  as 


u(t)  = 


m 

T  +  q 


(2.38) 


which  is  called  the  normalized  bubble  radius.  Initially  when  t  =  to. 


u(to)  =  1. 


(2.39) 


This  normalized  bubble  radius  can  be  applied  to  all  further  equations  involving  t  and  its 
derivatives.  Since  its  initial  value  is  unity,  it  makes  the  scaling  of  the  solutions  much  easier 
to  handle.  Therefore  the  normalized  bubble  radius  will  be  employed  in  those  results  which 
form  the  basic  ingredients  for  subsequent  analyses  of  bubble  dynamics. 


Another  important  characteristic  of  the  flaccid  bubble  problem  is  the  bubble  growth 
rate.  The  flaccid  bubble  growth  rate  is  defined  as  the  first  derivative  of  the  flaccid  bubble 
radius  with  respect  to  real  time,  t  This  rate  is  found  by  differentiating  the  cubic  flaccid 
bubble  equation.  Thus  from  (2.27) 

fw,[Cp(t)  +  K]«t)f  +  +  2f  =  0  .  (2.40) 

dC,  dr 

Substituting  for from  (2.17)  and  solving  for  one  finds  the  flaccid  bubble  growth 


rate  to  be 


d,  -  C„(.) 

*  '  1  +  ^r(t)  W,[Cp(t)  +  K] 


(2.41) 
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This  rate  is  a  function  of  the  flaccid  bubble  radius,  the  flow  parameters  -  Cp,  K,  Wr,  the 
free  stream  velocity,  v0,  and  also  some  body  parameters.  Particular  to  the  body  is  the  char¬ 
acteristic  diameter,  D,  and  its  pressure  coefficient  curve.  This  curve  leads  to  the  forcing 
function  on  to  a  given  body  that  ultimately  stimulates  the  bubble  nucleus  vaporous  growth. 
This  again  is  similar  to  the  result  achieved  by  Baker  [2]  and  as  he  shows,  it  relates  well  to 
experimental  results. 

To  normalize  this  result  with  the  same  parameter  used  in  the  flaccid  bubble 
radius  result,  equation  (2.38)  and  its  derivative  is  applied  to  (2.41)  resulting  in 


du 

dt 


,  Vn  , _ (  dC„^ 

fo(t)Wr^V(l  +q)(l  -  Cp(t 3TJ— gf 


1  +  |u(t)Wr>mrq[Cp(t)  +  K] 


(2.42) 


This  is  the  normalized  flaccid  bubble  growth  rate.  Its  initial  value  when  t  =  0  will  be 
discussed  in  a  later  section. 


2.3.4  Cakulatmofthe  Initial  Flaccid  BubbkGrowth  Rare,  ilQ) 

The  initial  growth  rate  begins  when  t  =  to  and  t,  =  x  =  0,  and  at  the  point  on  the 

body  where  Cp  =  -K.  Using  this  relationship  and  substituting  (2.37)  into  (2.41),  one 
gets  an  expression  for  the  initial  bubble  growth  rate  with  respect  to  real  time  which  is 


dr 

dt 


t  =  0 


-  i(l+q)Wrgvi  +  K 


(2.43) 


It  is  desired  to  analyze  the  growth  rate  in  terms  of  the  dimensionless  slow 
laboratory  time,  t,.  To  express  the  initial  growth  rate  in  terms  of  this  time  scale ,  the 


derivative  of  the  relationship 


(2.44) 


is  substituted  into  (2.43),  resulting  in 


dr 

dt. 


t,  ■  0 


(2.45) 
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Relative  to  this  time  scale  the  initial  bubble  growth  rate  loses  its  explicit  relationship  to  the 
characteristic  body  diameter  and  free  stream  velocity,  although  the  square  of  the  free  stream 
velocity  is  contained  in  the  Weber  number. 

Writing  the  initial  bubble  radius  with  respect  to  slow  time  in  normalized  form,  one 
needs  the  derivative  with  respect  to  slow  time  of  equation  (2.38).  Applying  this  to  the 
result  in  equation  (2.45),  one  sees  that  it  becomes 


du 

dts 


,0“|w,V(l+q)(l  ♦  K>(-d-^) 


Cp  =  -  K 


(2.46) 


This  is  the  normalized  initial  bubble  growth  rate  with  respect  to  the  slow  laboratory 
time.  One  gets  the  same  result  from  evaluating  equation  (2.42)  at  ts  =  0  and  convening  it 

to  the  slow  time  scale  via  equation  (2.44). 

So  far  the  two  initial  conditions  and  the  forcing  function  have  been  found  for  the 
Rayleigh-Plesset  equation  in  terms  of  the  two  dimensionless  time  scales  involved  in  the 
perturbation  expansion  for  the  point  when  Cp  =  -  K.  These  can  now  be  used  to  help  solve 

the  expanded  equations  which  will  be  derived  from  perturbation  theory. 


2.4  Formulation  of  the  Dynamical  Bubble  Equations  to  Order  e3 

The  pulse-forced  Rayleigh-Plesset  equation  (1.1)  written  in  terms  of  the  dimension¬ 
less  radius,  r,  with  respect  to  real  time,  t,  is 


(2.47) 


where  F(t)  is  some  forcing  function  characteristic  of  the  governing  body  and  the  flow 
properties1 .  As  noted  in  section  1.2  above,  this  is  the  governing  equation  for  vaporous 
growth  and  collapse  of  a  spherical  isothermal  cavitation  bubble.  It  is  a  second  order 
ordinary  differential  equation  requiring  die  two  initial  conditions  given  above  in  equations 


1  See  Appendix  D  for  the  derivation  of  this  dimensionless  form  of  equation  (1.1). 
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(2.37)  and  (2.45).  To  derive  an  appropriate  solution  to  this  equation  we  shall  expand  it  in 
ascending  powers  of  the  small  perturbation  parameter  e. 

The  method  of  multiple  scales  is  used  because  of  the  two  widely  varying  time  scales 
present  in  the  problem.  In  the  two-scale  form  [1 1],  the  small  parameter  e  allows  the  effect 

of  the  fast  bubble  time  to  be  pulled  into  the  slower  laboratory  time  in  a  uniform  fashion. 

The  slow  laboratory  time  scale,  ts,  is  the  characteristic  time  for  the  bubble  to  pass  along  the 


low  pressure  region  in  the  boundary  layer  of  a  particular  body.  This  time  scale  character¬ 
izes  the  forcing  function  since  it  depends  on  the  body.  The  fast  bubble  time  scale,  x,  is  the 
characteristic  time  of  the  free  stream  nucleus  oscillation  period.  This  scale  is  very  fast 
compared  to  the  slow  time  scale,  so  the  two  are  related  according  to, 

tS  =  et  .  (2.48) 

From  these  two  time  scales  and  the  small  perturbation  parameter,  the  dynamical  equation 
can  be  expanded  to  order  £3. 

We  shall  write  the  first  and  second  derivative  expansions  in  terms  of  the  two  time 
scales,  because  it  is  assumed  that 


r  =  r(ts,x,E). 


(2.49) 


Taking  the  derivative  of  r  with  respect  to  real  time,  using  the  chain  rule  and  (2.48),  one  has 


dr  _  dr  ^  \  dr 
dt  ~  dt,  +  e  dx  ’ 


(2.50) 


and  differentiating  again,  one  gets 

d2r  d2!  ,  2  d2r  (2.51) 

dt2  atj  edtgdx  e2  dx2 


Next  we  applying  the  above  first  and  second  derivative  expansions  to  a  general 
perturbation  expansion  of  r  of  the  form 

r  *  r0  +  rjE1  +  T2E2  +  r3E3  + ... ,  (2.52) 

which  leads  one  to  die  completely  expanded  final  forms.  They  are,  up  to  order  e3. 
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dr 

1  dro  ^ 

fro  <*r! 

dr,  dr2 

dr2  dr3 

dt  = 

e  dx 

dt*  +  dx 

+  £ 

K  +  d7 

+  e2 

dt,  +  dx 

(2.53) 


and 


d2r 

dt2 


1  ’ 

i  __ 

e2 

l 

QJ 

"~K> 

1 _ 

T 

e 

d2^ 


+ 


a2r0 

x 


+ 


a2!-,  a2r0 
dtsdT  +  at2 

s 


+  e 


+ 


a2r2  a^ 
atsdT  +  ats2 


(2.54) 


Determination  of  each  term  of  equation  (2.47)  can  be  based  on  the  general 
expansion  for  r  and  its  derivatives  (2.53)  and  (2.54).  As  was  shown  previously  in  section 
2.2,  the  forcing  function  F(et)  is  involved  only  in  the  first  order  of  e.  This  fact  leads  to  the 
first  order  equation  being  the  only  nonautonomous  expansion  equation. 

The  two  initial  conditions,  r(tg)  and  dr(to>/dt,  must  also  be  expanded  by  this 
method.  This  is  accomplished  by  equating  the  flaccid  bubble  initial  radius  of  (2.37)  and  the 
flaccid  bubble  initial  growth  rate  of  (2.45)  to  (2.52)  and  (2.53)  respectively,  and  equating 
like  powers  of  e. 

To  unify  the  non-dimensional  radius  we  introduce  from  equation  (2.38),  the 
substitution  of 


r  =  u  VT  +  q  ,  (2.55) 

where  q  =  KWr/4,  as  before.  This  comes  from  the  flaccid  bubble  radius  and  allows  one 

to  write  the  normalized  isothermal  Rayleigh-Plesset  equation  with  initial  conditions 
consistent  with  those  used  by  Baker  [2].  Once  solved  for  their  respective  uj's,  the  resulting 

set  of  differential  equations  and  initial  conditions  can,  using  (2.55),  be  substituted  into 
(2.52)  in  order  to  achieve  a  final  perturbation  solution  for  the  non-dimensional  bubble 
radius  as  a  function  of  real  time.  The  normalized  differential  equations  and  initial  conditions 
up  to  order  e3  are  displayed  below. 
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d2up  3  fdup  f _ ^ 

U°3t2  +  =  u’(l  + 


q)3/2  u0(l+q)3/2 


(2.56) 


with  initial  conditions. 


u0(0,0)=  1 


8u0(0,0) 


=  0 


(2.57) 


92ui  3_^q(M]  Ui  92u0 

dx2  +  up  (9x  J  +  u0  9x2 


up(l  +  q)3/2  uj 


=  -2 


92up  3  1  dup  dup  F(e,  t, 
axdts  +  2  up  dx  +  up (1  +  £ 


fs) _ 

q)3/2 


with  initial  conditions, 

M0.0)  =  0 

dui(O.O)  3u0(0,0)  1 


=  VWrV(l  +q)(l  +K) 


(-'-&) 


CP  = 


(2.58) 


(2.59) 
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e2_Qrder  Equation 
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(2.60) 


(2.61) 


.  (2.62) 


with  initial  conditions 


113(0,0)  =  0 

8u3(0,0)  9u2(0,0)  r  <2-63> 

dx  +  8ts  "  0 

2.5  Estimates  of  Flow  and  Perturbation  Parameter  Ranges 

In  the  previous  sections  many  parameters  have  been  defined.  It  is  of  prime  concern 
in  this  investigation  to  achieve  results  which  are  directly  related  to  the  characteristics  of  the 
fluid  flow  and  the  governing  body.  As  will  be  shown  later  all  flow  characteristics  are 
combined  into  one  parameter,  a  so  called  vortex  point  affix  defined  in  the  bubble  radius 
versus  radial  velocity  phase  plane  of  the  zero-order  solution.  This  parameter,  called  uv,  is 
used  throughout  the  analysis.  The  objective  of  this  section  is  to  determine  ranges  of  the 
parameters  that  define  the  flow  from  typical  values  of  obsevable  flow  quantities. 

There  are  three  basic  dimensionless  parameters  that  define  the  flow.  These  are  the 
radial  Weber  number,  Wr ,  the  cavitation  number,  K,  and  the  air  content  parameter,  y.  The 

radial  Weber  number,  as  defined  in  equation  (2.26),  is  based  on  the  initial  radius  of  the 
bubble  nucleus  in  the  free  stream.  For  water,  experiments  have  shown  that  typical  nucleus 
radii  vary  in  size  over  the  range 

5  <  Ro  <  150  nm  ,  (2.64) 

and  the  smaller  the  nucleii  the  greater  their  numbers.  Along  with  the  density  and  the  surface 
tension  for  water  at  70°  F,  the  Weber  number  also  contains  the  free  stream  velocity,  V0. 

Typical  velocities  obtainable  in  water  tunnels  vary  in  the  range 

0  <  V0  <  80  ft/s  .  (2.65) 

From  these  ranges,  we  find  the  range  estimate  for  the  radial  Weber  number  for  water  ro  be 
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The  cavitation  number,  as  defined  in  equation  (2.7),  depends  on  the  difference  in 
free  stream  static  pressure  and  vapor  pressure  as  well  as  the  speed  of  the  flow.  Previously, 
in  his  forced  bubble  growth  analysis,  Baker  [2]  used  a  small  range  of  cavitation  number 
values  typical  of  a  hemispherical  headform.  But  since  cavitation  numbers  can  reach  as  high 
as  2  or  3  for  certain  conditions,  for  this  estimate  we  use  an  expanded  range  that  could 
include  a  variety  of  models.  This  range  varies  from 

0.1  <K  <3.0.  (2.67) 

The  Weber  and  cavitation  numbers  are  both  contained  in  the  parameter  q,  as  defined 
in  equation  (2.32).  Using  the  ranges  for  these  in  equations  (2.66)  and  (2.67),  we  find  the 
range  of  q  to  be  approximately 

0  <  q  <  911  .  (2.68) 

The  air  content  parameter  depends  on  the  size  of  the  initial  nucleus  radius  and  the 
partial  pressure  of  the  gas  at  the  liquid-gas  interface,  pa.  This  partial  pressure  is  found 

using  Henry's  Law  where 


Pa  =  •  (2.69) 

Henry's  Law  constant,  p,  is  in  general  a  function  of  pressure  and  temperature.  But  for 
pressures  well  above  atmospheric  it  is  nearly  independent  of  pressure,  but  it  still 
depends  on  the  temperature.  At  a  temperature  of  70°  F,  Henry's  Law  constant  is 

»  -  °"  £=  •  (2.70) 

where  “ppm”  refers  to  molal  parts  per  million.  The  amount  of  any  dissolved  gas,  a,  in 
water  tunnel  tests  typically  varies  within  the  range 

3  <  a  <  15  ppm  .  (2.71) 

The  air  content  parameter  is  defined  as  the  ratio  of  the  saturation  air  pressure  to 
the  pressure  due  to  surface  tension,  namely: 
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Y  = 


paR0  paR0 


(2.72) 


2o  2o 

Substituting  the  range  of  variables  from  (2.64),  (2.70)  and  (2.71)  into  (2.72),  one  finds  the 
range  for  the  dimensionless  air  content  parameter  to  be  approximately 


1  <  Y  <  105  .  (2.73) 

The  relationship  between  the  slow  laboratory  time,  tj,  and  the  fast  bubble  time,  t,  is  given 

by  the  small  perturbation  parameter,  e.  It  can  be  expressed  as 


R0  o 

£  =  irw?’ 


(2.74) 


which  directly  combines  the  key  flow  parameters  in  the  Weber  number  and  the  characteristic 
body  diameter,  D,  into  the  single  parameter,  e.  Figure  2.6  below,  shows  the  dimensionless 
relationships  between  the  perturbation  parameter,  e,  the  body,  the  nucleus  sizes  and  the 
measurable  flow,  and  fluid  properties  in  accordance  with  equation  (2.74). 
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The  ratio  of  "bubble  time"  to  "laboratory  time"  equals  e. 


0.001  0.002  0.004  0.007 


Figure  2.6  Values  of  Perturbation  Parameter,  e,  for  Ranges  of  Radial  Weber 
Number  and  Nucleus  Size  to  Characteristic  Body  Length. 
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CHAPTER  3 

THE  ZERO-ORDER  SOLUTION 

3.1  The  First  Integral 

From  the  previous  development  the  zero-order  form  of  the  Rayleigh-Plesset 
differential  equation  in  terms  of  the  normalized  dimensionless  bubble  radius  u0(t)  is 


32u0  3  fdu0 

U°9x2  +  2  dt 


3  ("duo?  _ X 

2Vd<cJ  u3(l  + 


(1  +  q)5/2  u0(l+q)3/2 


with  the  initial  conditions 

uo(0)  =  1  (3.2) 

dun 

af  t.0  =  0  '  (33) 

It  is  assumed  at  this  point  that  the  zero-order  equation  is  independent  of  the  slow  laboratory 
time,  tj,  and  therefore  can  be  written  as  an  ordinary  differential  equation  in  terms  of  the  fast 

bubble  time,  x.  The  reasons  for  this  independence  are  addressed  in  section  3.5  below. 

An  important  fact  about  this  zero-order  equation  is  that  it  is  autonomous  and 
therefore  its  first  integral  is  easily  found.  By  letting 


v  "  dx  ’  VA> 

one  can  write  (3.1)  as  two  coupled  first  order  equations,  namely  equation  (3.4)  and 

d(uov2)  =  2y  2  u0  (3.5) 

duo  u0(l+q)5/2  (1  +  q)3/2 

In  anticipation  of  the  discussion  in  section  3.5  we  integrate  both  sides  of  equation  (3.5) 
partially  with  respect  to  x  in  order  to  obtain  a  first  integral,  one  finds  it  to  be 

U”V*  =  (1  +2q)3'jl"<Uo>  '  (1  +  q)3'2  +  [A<t,)  +  (1  +  q)5«]  ’  <3'6> 
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where  the  quantity  ( ts)  +  is  the  constant  of  integration  which  must  satisfy 

the  initial  condition  of  equation  (3.2).  To  also  satisfy  the  condition  in  (3.3)  at  ts  =  t  =  0, 
however,  A(0)  must  equal  zero,  which  leaves 


UqV2  = 


1 


(1  +q) 


3/2 


(TM)ln(Uo) 


Uq  +  1 


(3.7) 


Further  arguments  that  A  =  constant  =  0  and  that  the  zero-order  equation  can  be  indepen¬ 
dent  of  ts  will  be  shown  in  section  3.5. 

The  left  hand  side  of  this  first  integral  is  proportional  to  the  kinetic  energy  of  the 
cavitation  bubble  motion,  whereas  the  right  hand  side  is  proportional  to  the  potential 
energy,  -  V,  of  the  bubble.  Thus  equation  (3.7)  can  be  plotted  in  a  potential  energy  versus 
Uq  plot  and  in  a  phase  plane  of  &o  versus  Uq  as  shown  in  Figure  3.1  and  the  nature  of  its 
phase-plane  trajectories  analyzed. 


3.2  Phase-Plane  Analysis 

Had  the  quantity  A  not  been  zero  in  equation  (3.6),  suggesting  a  step-forced  zeroth 
order  equation,  the  plot  of  the  potential  function  and  its  corresponding  phase  plane  trajecto¬ 
ries  would  lead  to  the  well-known  solution  of  the  autonomous  equation  arrived  at  by  Baker 
[2]  and  others  before  him.  The  phase  plane  then  includes  both  vortex  and  saddle  points, 
and  a  separatrix  which  separates  closed  periodic  trajectories  from  ones  that  grow  without 
bound. 

Under  the  assumption  that  A  is  equal  to  zero,  equation  (3.7)  can  be  plotted  in  the 
phase  plane  and  its  trajectories  analyzed.  Figure  3.1  shows  schematically  the  results  of  this 
analysis  for  a  typical  value  of  the  dissolved  air  content  parameter,  y  =1.4,  after  Baker  [2]. 
The  saddle  point  vanishes  and  all  trajectories  originate  from  the  point  at  Uq  =  1 ,  v  =  0.  In 

figure  (3.1)  the  expression  (1  +  Q)  is  simply  an  approximate  form  of  the  exact  expression 
for  Vl  +  q  used  in  this  analysis.  Baker  assumed  that  q  was  a  small  number  less  than  one 
which  is  not  necessarily  the  case  in  this  investigation.  Despite  this  difference,  figure  3.1 


NORMALIZED  BUBBLE  WALL 
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DIMENSIONLESS  INITIAL  RADIUS,  r (0)  =  1  +  Q  (Q<  1) 


Figure  3. 1  Schematic  Diagram  of  the  Phase  Plane  and  Level  lines  of  Potential 

Energy  Plot  for  the  Zero-Order  Analysis,  after  Baker  [1]. 
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shows  the  same  qualitative  phase  plane  relationships  as  found  here.  It  is  seen  that  the 
location  of  the  minimum  potential  energy  point  for  a  given  value  of  q  corresponds  to  the 
vortex  point  in  the  phase  plane  about  which  the  trajectories  are  focused.  These  vortex 
points  move  from  left  to  right  as  the  value  of  q  decreases  to  zero.  Free  oscillations  of  the 
bubble  will  then  have  amplitudes  determined  by  distance  of  the  minimum  energy  level 
below  zero  because  the  system  being  conservative,  the  initial  energy  is  the  total  available 
and  it  is  set  by  the  initial  conditions  to  be  at  zero.  The  trajectories  on  the  right  represent 
oscillations  of  expansion  with  significant  vaporous  growth  accompanied  by  a  reduction  in 
the  partial  air  pressure  in  the  bubble.  Those  on  the  left  are  compressive  oscillations  in 
which  the  air  contained  in  the  bubble  supplies  the  restoring  pressure  as  before  and  the  mass 
of  vapor  in  the  bubble  is  reduced. 

The  dotted  curve  in  the  upper  potential  energy  plot  corresponds  to  the  case  where 
the  trajectory  and  the  vortex  are  isolated  at  the  origin  in  the  phase  plane.  The  motion  in  this 
critical  case  is  null  except  in  the  possible  event  of  a  momentary  random  disturbance  which 
would  excite  free  oscillations  of  the  bubble.  Trajectories  for  this  case  would  be  centered 
about  the  vortex  point  at  Uo(0)  =  1,  v(0)  =  0.  The  oscillation  amplitude  would  depend  on 
the  intensity  of  the  momentary  pulse.  It  would  supply  an  initial  increment  of  energy  to  the 
bubble  causing  the  initial  energy  level  in  figure  3.1  to  move  from  zero,  as  required  by  the 
present  initial  conditions,  to  some  higher  energy  level.  These  trajectories  are  not  shown  in 
figure  3.1  because  a  random  disturbance  for  that  particular  case  is  not  contemplated.  We 
require  the  initial  conditions,  Uo(0)  =  1,  v(0)  =  0,  to  be  satisfied  strictly.  Nevertheless  the 
natural  frequency  of  these  vibrations  is  of  interest  and  these  are  discussed  for  small 
oscillations  in  the  section  3.41  below. 

The  importance  of  this  critical  case  is  that  it  separates  the  bubble  modem  into  the 
distinct  types,  oscillations  of  compression  and  expansion.  Trajectories  to  the  right  of  the 
critical  point  represent  larger  amplitude  oscillations  involving  vaporous  growth  and  a  small 
amount  of  air  induced  motion.  Since  this  investigation  considers  primarily  the  vaporous 
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growth  region,  the  trajectories  to  the  right  of  the  critical  point  will  be  emphasized,  although 
as  will  be  found  later  there  can  be  some  flow  conditions  under  which  those  on  the  left  can 
play  a  role. 

The  location  of  the  vortex  points  corresponding  to  the  minimum  of  the  potential 
energy  curves  is  found  by  setting  the  first  derivative  of  equation  (3.7)  equal  to  zero  and 
solving  for  Uq,  resulting  in  the  affix  of  the  vortex  point  being  defined  analytically  by 

u.  =  Vdhi  •  (3-8> 

Then  for  any  air  content,  at  the  critical  point  where  Uy  =  1,  a  critical  value  of  q  can  be 

defined  as 


qCrit  =  Y  -  1  •  (3.9) 

If  q^j  =  0,  then  y  =  1  which  sets  a  minimum  value  for  the  dissolved  air  content  parame¬ 
ter  for  positive  q  values.  Since  q  depends  on  the  cavitation  number  and  the  Weber  number, 
a  limiting  value  of  q  =  0  corresponds  to  K  =  Oorv0  =  0.  This  result  is  equivalent  to  that 

achieved  by  Baker  [2],  except  in  this  case  the  value  of  q  is  not  necessarily  less  than  one. 
Therefore,  the  denominator  of  equation  (3.8)  cannot  be  approximated  using  the  binomial 
expansion  as  stated  before. 

3.2.1  Measurable  Hfiffi  Parameters  and  the  Value  sL  uv 

The  analytic  expression  for  the  vortex  point  affix,  Uy,  in  equation  (3.8)  contains  all 
three  of  the  dimensionless  flow  parameters  described  in  section  2.5,  namely:  y,  Wp  and 
K,  but  it  does  not  include  the  body  size,  D,  because  it  applies  only  to  the  bubble.  By  using 
Uy  and  the  small  perturbation  parameter,  e  which  does  not  include  the  air  content  parameter, 
y,  or  the  cavitation  number,  K,  but  all  flow  variables  and  headform  sizes  can  be  contained 
in  e  and  Uy.  Therefore  it  is  useful  to  complement  the  relationships  of  figure  2.6  with 
another  illustration  that  includes  the  remaining  physical  variables.  We  shall  consider  cases 
of  Uy  2, 1  only. 
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The  range  of  cavitation  numbers  is  constrained  by  likely  values  of  Cplmjn  for 
various  submerged  bodies.  Thus  if  a  particular  vortex  point  is  desired,  various 
combinations  of  flow  variables  and  headfonn  sizes  can  be  used  to  arrive  at  this  number. 
Using  the  critical  values  of  q  and  the  air  content  parameter  arrived  at  above,  along  with  the 
ranges  of  values  for  the  Weber  Number  and  cavitation  number  from  section  2.5,  one  sees 
that  values  of  uv  vary  in  the  range  from 


1  ^  My*  10  •  (3.10) 

A  relationship  for  these  parameters  is  shown  in  figure  3.2.  The  plot  parameter,  z,  in  the 
figure  is  used  to  relate  the  two  scales  and  is  given  by 


1  +  q  >/y 


(3.11) 


To  use  the  graph,  a  cavitation  number  is  chosen  and  traced  horizontally  to  the  appropriate 
Weber  number.  This  point  corresponds  to  a  certain  value  of  z  directly  below  it.  Finding  the 
vertical  intersection  of  this  value  of  z  with  the  chosen  air  content,  y,  one  can  then  move 


horizontally  to  the  scale  on  the  right  in  order  to  find  the  affix  of  the  vortex  point .  Finally, 
since  the  radial  Weber  number  is  known  one  only  needs  to  know  the  ratio  Rq/D  from 
figure  2.6  in  order  to  find  the  proper  value  of  £  for  the  hydrodynamic  configuration.  And 
since  uv  is  now  known,  values  of  u^m  and  X  follow  from  figures  3.3  and  3.4  below. 


The  phase  plane  analysis  shows  that  the  first  integral  yields  closed  trajecto¬ 
ries.  Therefore  both  the  zero-order  radius  and  growth  rate  are  periodic  with  respect  to  x. 
By  separating  the  differential  dx  from  Uq  and  its  differential  in  equation  (3.7),  one  has 


35  30|25  2 
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Because  the  solution  is  periodic  this  integral  can  be  evaluated  over  half  a  period,  say  T/2, 
and  its  evaluation  takes  the  form, 


(3.13) 


1 

where  the  upper  limit  on  the  right  hand  integral,  um(uv),  is  the  maximum  normalized  radius 
the  oscillating  bubble  will  reach. 

The  bubble  grows  from  it's  normalized  value  of  unity  to  this  maximum  value  and 
then  returns  to  unity.  The  maximum  occurs  at  exactly  half  the  period  of  oscillation  and 
corresponds  to  half  way  around  a  trajectory  in  the  phase  plane  plane.  The  value  of  u„,  as  a 
function  of  uv  is  obtained  from  solving  equation  (3.7)  for  u0  >  1  under  the  condition  that 


v  =  0,  thus 

2u.21"(0  •  <&  *  1  =  0  . 


(3.14) 


We  see  by  inspection  when  u<)  =  um  =  1  for  any  uv ,  that  equation  (3.14)  is 
satisfied.  On  the  other  hand  if  uv  =  1  then  um  =  1  is  the  physically  acceptable  root  of 
equation  (3.14).  Therefore  at  the  initial  point  where  uq(0)  =  1  there  will  be  no  bubble 
motion  if  uv  =  1,  as  discused  in  section  3.2  above.  Although  a  closed  form  of  the  required 
solution  for  um(uv)  is  unobtainable,  equation  (3.14)  can  easily  be  solved  for  its  inverse, 
u^u,,,).  But  here  we  have  used  the  numerical  methods  of  regula  falsi  or  bisection  in  order 
to  get  the  desired  results.  Figure  3.3  shows  u^  as  a  function  of  uv  and  includes  a  nonlinear 
least  squares  curve  fit  over  the  range  of  u*  values.  The  curve  fit  matches  the  calculated  data 
very  well  but  it  was  not  used  in  the  numerical  calculations  in  order  to  secure  the  greatest 
accuracy  in  subsequent  numerical  work  in  which  values  of  uv  must  be  prescribed. 

Defining  the  integral  of  equation  (3. 1 3)  in  terms  of  a  prescribed  uv  from 


max 
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Curve  fit  for  umax. 


Figure  3.3  Maximum  Normalized  Zero-Order  Bubble  Radius  for  a  Range  of  uv. 
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equation  (3.14)  and  excluding  the  constant  term  in  front  of  the  integral,  we  write  the 
half-period  integral  as 


“m 


IK)  = 


*\j 2uJ  In 


IL 


(O  -  C2  +  i 


dC  • 


(3.15) 


Since  the  half-period  integral  is  now  known,  we  can  write  the  oscillation  period  of  the 
solution  u0(uv,t)  to  be 


T(uv)  =  2(1  +  q)3/4I(uv)  .  (3.16) 

The  period  is  a  function  of  both  uv  and  the  flow  parameters  contained  in  the  quantity  q. 

By  using  this  period  the  non-dimensional  bubble  time,  x,  can  also  be  normalized. 
This  is  done  by  defining  the  normalized  dimensionless  bubble  time,  x,  as 

X 

x  =  f  *  (3-17) 

When  x  =  0,  x  =  0  and  after  one  period  when  x  =  T,  x  =  1.  This  also  allows  the  value 
of  Um  to  occur  at  x  =  1/2. 

3-4  Zero-Order  Equation  Solution 

From  the  first  integral  of  equation,  (3.7)  and  the  normalized  bubble  time  from 
equation  (3.17),  the  exact  first  integral  of  the  zero-order  equation  now  looks  like 

.  / 2ujln(u0)  -  uj  +  1 

duo  _T__y\  /  _ _  (3.18) 

^‘(l+q  )3M\  4 

If  one  separates  variables  in  equation  (3.18)  the  problem  of  finding  the  inverse  solution, 
x(uv,Uo),  is  reduced  to  the  quadrature, 
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(3.19) 


where  the  period  parameter ,  X,  is  a  constant  for  a  given  value  of  uv  and  is  related  to  the 
period  by 


,,  ,  T(Uv)  _  , 

X  v  =  (1  +q)3/4  =  2I(Uv)  •  (3.20) 

These  results  make  the  zero-order  solution  a  one-parameter  family,  u0(x,uv),  which  is 

convenient  for  further  analysis.  It  is  noted  that  the  other  major  components  involved  in  this 
equation  also  have  been  derived  with  this  dependent  parameter,  such  as  um(uv)  and  X(uv). 
Figure  3.4  shows  a  graph  of  the  period  parameter  over  the  range  of  uv  along  with  a  least 
squares  curve  fit  of  the  function.  Next  we  shall  adopt  the  notation  of  equation  (3.20)  and 
use  it  in  equation  (3.18)  in  order  that  it  will  have  the  same  notation  as  equation  (3.19) : 


w  -  *“v> 


2ln(u0)  -  uj  +  1 


(3.18  a) 


Equations  (3.4)  and  (3.18a)  have  been  used  in  order  to  calculate  a  number  of  phase 
plane  trajectories  for  the  range,  0.50  <1  uv  £  1.50,  enabling  us  to  see  how  the  trajectories 
change  as  the  value  of  uv  passes  through  unity.  These  trajectories  are  shown  in  figure  3.5 


below. 

A  solution  to  equation  (3.19)  would  yield  x  =  x(uo).  By  approximating  the  logarithmic  air 

content  parameter  with  a  cubic  polynomial,  Baker  [2]  found  an  approximate  solution  for 
x(uq)  in  terms  of  elliptic  integrals  and  functions.  Of  more  general  use  is  the  inverse  of  this 
function,  or  uq  =  uq(x).  We  shall  proceed  numerically  because  Baker’s  solution  holds  for 
a  rather  small  range  of  uv  and  it  can  not  be  inverted 


3.4.1  Zero-Order  Small  Oscillation  Study 

From  the  phase  plane  plot  of  u©  vs  uq  it  is  seen  that  for  a  given  value  of  uv  that  the 
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Curve-fit  and  integrated  values  of  X  compared. 


Figure  3.4  Values  of  the  Period  Parameter,  X,  for  a  Range  of  uv  Values. 


44 


0.0  0.5  1.0  1.5  2.0 


Normalized  radius,  u(x). 


Figure  3.5  Phase  Plane  Trajectories  Calculated  for  Several  uv 
Values  in  the  Range  0.5  is  uv  ^  1.5. 
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curves  are  closed  and  nearly  elliptical,  especially  for  small  values  of  Uy  near  1.0.  Conse¬ 
quently,  an  ad  hoc  small-oscillation  form  or  the  zero-order  solution  was  sought  for  the  case 
when 

uv  =  1  +  8.  (3.21) 

where  8  «  1.  For  this  case  an  approximate  form  of 
the  first  integral  was  sought  using  the  equation  for  an 
ellipse.  To  most  closely  approximate  the  form  of  the 
first  integral  the  ellipse  was  chosen  to  be  centered  at  c, 
one  -half  the  distance  between  the  two  points,  1  and 
um,  on  the  u0  axis.  The  length  of  the  vertical.axis  is 
determined  by  the  requirement  that  the  period  of 
oscillation  must  also  agree  with  the  correct  result.  The 
use  of  these  conditions  leads  to  an  equation  for  an  ellipse  of  the  form 

.2 

U°  fu0  -  c¥  .  (3.22) 

b2  +  l  *  J  =  ' 

where 

a  =  \  (“m  -  !)  *  (3.23) 

b  =  (uo)max  .  (3.24) 

and 

c  =  2  +  •  (3.25) 

The  coefficients  a,  b  and  c  are  all  functions  of  uv  since  u,n  =  um(uv)  and  b  comes  from  the 
first  integral  equation  which  is  also  a  function  of  uv.  The  quantity  b  equals  the  maximum 
value  of  Uq,  which  must  be  arrived  at  numerically.  But  as  will  be  shown  below,  the 
approximate  solution  can  be  written  independently  of  the  parameter  b. 
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When  this  result  is  used  in  equation  (3.28)  and  the  normalized  bubble  time  from 
equation  (3.17)  is  applied,  the  small  oscillation  approximate  solution  becomes 

Uo(Uv,x)  =  c  -  a  cos(2tc  x)  .  (3.31) 

Equations  (3.23)  and  (3.25)  show  that  c  and  a  depend  on  um(uv),  and  so  equation  (3.31) 
depends  on  u* ,  as  well  as  the  normalized  bubble  time,  x.  As  with  the  exact  solution,  the 

approximate  solution  has  no  explicit  dependence  on  the  parameter  b.  Equation  (3.31)  is 
designed  to  have  the  correct  nonlinear  amplitude  and  period  as  determined  by  Uy,  however. 
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Moreover,  when  x  =  0  equation  (3.31)  yields  Uo(Uy,0)  =  1.  Differentiating  (3.31)  with 

dlln 

respect  to  x  and  evaluating  it  at  x  =  0  one  has 


x  =  0 


=  0.  Thus  equation  (3.31) 


satisfies  the  initial  conditions. 

Next  we  need  to  find  the  variation  of  u„,  with  uv,  when  5  =.uv  - 1  «  1  from 
equation  (3.21).  Knowing  that  also  doesn't  vary  greatly  from  unity  for  small 
oscillations,  we  can  be  approximate  it  by 


um  =  1  +  ym  .  (3.32) 

where  ym  «  1,  but  ym  is  slightly  larger  than  6.  First  one  substitutes  equation  (3.32)  into 

a 

equation  (3.14)  and  expands  the  logarithmic  term  to  0(ym)  in  order  to  obtain  a  quadratic 

equation  in  ym. which  can  be  solved  for  ym(uv).  Then  equation  (3.21)  can  replace  uv  in 
this  solution  in  order  to  expand  it  to  CK82).  Thus  one  gets  the  rather  accurate  expansion. 


u™  =  1  +  25  +  |82  . 


(3.33) 


Substituting  this  result  into  c  and  a  for  um  in  order  to  get  a  small  oscillation 
expansion  for  ug(x;  uv)  we  see  that  it  looks  like 


8 

u0  =  1  +  +  g)(l  -  cos2jcx)  .  (3  34) 

It  is  seen  that  as  8  -*  0,  both  Uq,  and  Um  -•  1  which  agrees  with  the  initial  condition. 

In  order  to  get  an  independent  evaluation  we  now  find  a  limiting  small-oscillation 
form  of  the  period  parameter,  X,  as  given  in  detail  in  Appendix  C.  Here  we  only  outline 
the  analysis.  To  this  end  we  may  generalize  equation  (3.32)  by  writing  uq  =  1  +  £  so  that 
when  £  =  ym,  uo  =  um,  but  generally  0  £  £  £  ym  because  1  <  uq  £  um  .  Then  the  dummy 
of  integration  in  equation  (3.15)  is  replaced  by  uq  and  the  denominator  inside  the  radical  of 
equation  (3. 15)  can  be  expanded  to  0(C3)  and  expressed  in  factored  form  using  the  three 
roots  of  the  approximating  cubic  at  the  three  £  values  of  0,  £  ym,  £  bo-  The  period  integral 
then  becomes 
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I(uv)  = 


fuv2C(ym-C)(b0*O 


(3.35) 


Additional  expansions  of  the  integrand,  as  explained  in  Appendix  C,  and  term-by-term 
integrations  of  the  result  give  the  respective  components  of  I  in  terms  of  ym .  The  sum  of 


these  components  gives  the  rough  estimate, 


(3.36) 


Because  of  equation  (3.35),  equation  (3.36)  contains  the  limiting  value,  b0  =  3  as  ym-»  0 

and  uv->  1.  Equation  (3.36)  shows  that  as  uv-»  1, 1  =-p ,  exactly.  But  since  b0  and  ym 

\2 

depend  on  uv  =  1  +  5,  one  can  replace  the  estimate  of  (3.36)  as  an  expansion  of  the  half¬ 
period  integral  in  powers  of  5  by  consistent  use  of  expansions  of  the  formulas  of  appendix 
C  for  ym  and  b0.  But  we  prefer  to  use  the  expansion  of  equation  (3.33).  This  expansion 
can  be  used  to  expand  the  quantities  ym  and  b0,  keeping  terms  up  to  order  S2,  in  order  to 

obtain  a  small  oscillation  form  of  the  half-period  integral,  replacing  equation  (3.36). 
Multiplying  the  result  by  a  factor  of  2,  the  small  oscillation  form  of  the  period  parameter  is 


X  =  W2  1  +  1 5  + 


(3.37) 


As  5  -»  0,  the  value  of  X  -*  JtV2  which  sets  the  value  of  X  for  uv  =  1.  As  explained  in 


section  3.2,  this  value  of  1  is  the  period  parameter  for  free  oscillations  which  are  expanding 
in  (me  part  of  the  cycle  and  are  compressive  during  the  other  part  It  is  also  noted  in  ap¬ 
pendix  C  that  equation  (3.35)  is  a  complete  elliptic  integral  of  the  third  kind,  from  which  the 
limiting  value  of  has  been  found  in  order  to  check  the  limit  from  equation  (3.36). 


These  small  oscillation  results  are  useful  for  similar  studies  of  the  first-order  equation.  A 
comparison  the  small  oscillation  results  with  numerical  calculations  shows  that  the  small 
oscillation  equations  are  very  accurately  valid  for  the  range  of  Uy  from  unity  up  to  1.03. 
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3.4.2  Numerical  Solution 

Since  there  is  no  closed  form  solution  for  uq(x),  an  integral  such  as  the  one  in 

equation  (3.19)  is  generally  evaluated  numerically.  For  this  task  we  have  used  Simpson's 
rule  for  incremental  steps  of  u^.  The  integrand  on  the  right  side  of  equation  (3.19)  contains 
integrable  singularities  at  Uq  =  1  and  um,  and  their  presence  must  be  accounted  for  in  the 
numerical  work. 

This  integrand  from  the  right  side  of  equation  (3.19),  labeled  F,  is 

(3.38) 

The  denominator  of  F  has  roots,  as  seen  from  equation  (3.14),  at  u0  =  1  and  Uq  =  um. 

By  using  a  small  perturbation  analysis  about  these  points  approximate  formulas  can  be 
derived  to  avoid  discontinuities. 

This  is  done  by  first  perturbing  Uq  about  the  point  Uq(0)  =  1,  so  we  let 
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Here  Fj  and  F2  are  equation  (3.38)  evaluated  at  x  equal  h  and  2h  respectively.  When  these 
two  equations  for  A  and  B  are  solved,  a  so  called  "start"  value  for  equation  (3.38) 
integrated  over  2h  is 

2yf2  r- 

start  =  3xaoh(4Fi  *  ^p2)  •  0.42) 

The  constant  X(Uy)  is  needed  to  make  the  'start'  value  consistent  with  equation  (3.19). 

This  is  labelled  a  'start'  value  since  it  is  used  at  the  beginning  of  the  integration  when 
x  =  0. 

The  other  singularity  occurs  at  the  point  where  Uo(^)  =  tv  Applying  the  same 
method  for 

uo  =  Um  -  5  ,  (3.43) 

one  finds  that  the  first  order  approximation  of  equation  (3.38)  reduces  to 

fUl-aV8  +  bs3«  .  (344) 

Solving  for  the  A  and  B  as  before  yields  a  so  called  "stop"  value  for  equation  (3.38)  over 
2hof 

St°P  =  153^jh(4V5Fl  '  Fz)  *  (3-45) 

Here  Fj  and  F2  are  equation  (3.38)  evaluated  at  x  equal  ^  -  h  and  |  -  2h,  respectively. 

This  is  labelled  a  "stop"  value  since  it  is  used  at  die  end  of  the  first  half  period  of  the 
solution  when  x  » 

The  fact  that  the  zero-order  solution  is  periodic  allows  the  "start"  value  to  be  used 

over  the  first  two  steps  of  the  integration  nearest  x  =  0  and  also  over  the  last  two  steps  at 

the  end  nearest  x  =  1.  Similarly,  the  "stop"  value  can  be  used  over  the  two  steps  just 
preceding  and  just  following  x  =  |.  Periodicity  also  warrants  integration  only  up  to  the 

x  =  ^  point  since  the  values  of  the  solution  between  one-half  and  one  will  correspond  to 
die  reverse  of  those  between  one  and  one-half. 
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The  numerical  calculation  of  the  zero-order  solution  involves  first  finding  a  value  of 
um  for  a  given  value  of  uv  using  equation  (3.14).  Simpson's  rule  is  then  employed  over 
the  range  of  steps  between  the  first  two  and  last  two  for  0  <  x  <  | .  The  curve  over  this 

interval  is  well  behaved  and  resembles  a  versine  curve,  as  exemplified  in  the  small 
oscillation  analysis.  Then  the  start  and  stop  values  are  added  to  the  result  These  values 
are  then  used  in  a  cubic  interpolation  scheme  to  produce  evenly  spaced  values  of  u0  over 

the  whole  range  of  x  from  0  to  1.  A  Fortran  source  code  for  these  calculations  is  given  in 
appendix  A. 

The  numerical  solution  for  u0  and  its  derivative  versus  x  for  various  values  of  uv 
are  shown  in  figure  3.6.  It  is  seen  that  the  curves  resemble  versine  curves,  insuring  a 
smooth  transition  to  the  small  oscillation  result.  To  get  a  better  picture  of  the  zero-order 
solution's  periodicity,  figure  3.7  shows  Uq  and  its  derivative  over  two  periods  in  x  for  a 
single  value  of  =  1.05  where  they  have  nearly  trigonometric  shapes.  These  results 
pertain  to  the  case  in  which  uv  >  1.  In  order  to  compare  such  solutions  with  solutions  when 
uv  <  1,  figure  3.8  shows  calculated  results  for  uv  =  1.5, 1.0  and  0.5.  These  curves  appear 
to  show  a  180°  phase  shift  as  uv  passes  through  unity.  The  cause  of  this  can  be  seen  in  the 
phase  plane  of  figure  3.5  in  which  trajectories  on  the  right  of  Uq  =  1  are  traced  in  the 
counter  clockwise  sense  as  x  increases,  starting  on  the  upper  part  of  a  loop.  Trajectories  to 
the  left  start  on  the  lower  branch  of  a  loop  and  also  proceed  in  the  counter  clockwise 
direction  with  increasing  x.  The  elongated  shape  of  the  trajectory  for  uv  =  1.5,  showing 
rather  modest  values  of  u  especially  on  its  right,  insures  that  in  figure  (3.8)  the  central  pan 
of  the  solution  near  x  =  1/2  will  be  relatively  gendy  rounded  compered  to  its  behavior  in  the 
neighborhood  of  Uq  =  1.  In  this  latter  region  figure  3.5  shows  larger  velocities  which 
increase  rapidly  from  zero  near  Uq  =  1.  Therefore  in  figure  3.8  the  solution  is  relatively 
steep  near  x  =  0  and  1.  On  the  trajectory  for  uv  =  0.5  of  figure  3.5  the  larger  and  rapidly 
changing  velocities  are  found  near  its  left  extremity.  Consequently,  figure  3.8  shows  sharp 
minima  at  x  *  0.5  and  1.5,  but  gende  rounding  near  its  maxima  at  x  =  0  and  x  =  1. 


Zero-order  solution  for  several  vortex  point  affixes,  u, 


Figure  3.6  Numerical  Zero-Order  Solution  as  a  Function  of 
Normalized  Bubble  Time  for  Several  Values  of  uv. 


0.0 


0.5  1.0 

Normalized  bubble  time,  x. 


1.5 


Figure  3.8 


Calculated  Oscillations  of  Expansion,  uv  >  1,  Compared 
with  Compressive  Oscillations,  uv  <  1. 


55 


3.4.3  Approximate  Solution  Using  Fourier  Series 

It  is  seen  that  the  small  oscillation  approximate  solution  resembles  the  first  two 
terms  of  an  even  function  expressed  in  a  Fourier  cosine  series.  This  encouraged  the 
development  of  a  compatible  approximate  solution  for  the  zero-order  solution  using  Fourier 
series  expansions. 

Considering  the  fact  that  the  zero-order  solution  is  periodic  and  also  an  even 
function  as  displayed  by  the  numerical  results,  it  seems  fitting  that  the  data  from  the 
numerical  solution  over  one  period  could  easily  be  fitted  to  a  truncated  Fourier  cosine  series 
in  which  the  a„  vary  continuously  with  uv.  Then  by  doing  this  for  a  large  number  of  uv 
values,  curve  fits  for  the  needed  a„(uv)  terms  can  be  found  using  a  least  squares  method. 

The  array  of  equations  for  the  appropriate  number  of  Fourier  cosine  terms  represents 
essentially  a  closed  form  approximate  Uq  solution  as  a  function  of  x  for  a  given  value  of  uv. 

The  advantages  of  having  such  a  solution  are  found  in  the  relative  speed  and  ease  of 
computation,  along  with  the  fact  that  many  of  the  complications  associated  with  numerical 
methods  are  avoided. 

The  first  thing  is  to  determine  the  number  of  terms  to  include  in  the  truncated  series. 
Initially  100  a„  coefficients  were  calculated  by  recursion  [14]  for  a  number  of  uv  values. 
The  same  was  done  for  50  a„  coefficients  and  the  results  compared  to  the  corresponding 
numerical  data.  It  was  seen  that  a  truncation  to  50  coefficients  loses  very  little  accuracy 
compared  to  the  numerical  results. 

Adopting  a  50  coefficients  truncation  for  the  Fourier  series  at  numerous  values  of 
Uy  over  the  range  from  1  to  10, 50  separate  curve  fits  were  then  applied  to  the  a^s  over 
their  corresponding  range  of  Uy.  This  is  to  say  that  only  the  lower  numbered  terms  were 
used  for  smaller  values  of  Uy,  and  as  Uy  increased  so  did  the  number  of  terms.  It  was 
found  that  all  except  a 0  were  effectively  handled  if  they  were  represented  by  a  rational 
fraction  in  terms  of  Uy.  A  quartic  polynomial  seemed  most  appropriate  for  a^  Each  curve 
fit  was  determined  using  a  least  squares  method.  Then  given  a  value  of  Uy  the  appropriate 
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number  and  value  of  a^s  can  be  determined  and  each  multiplied  by  it's  cosine  of  its 
frequency  and  summed  up  over  a  period  in  x  from  0  to  1,  yielding  a  truncated  series  zero- 
order  solution. 

The  convenience  of  curve  fitting  does  cost  some  accuracy.  It  was  found  that  in 
summing  up  the  cosine  terms  the  required  value  of  uo(0,uv)  =  1  was  missed  by  various 
small  amounts  depending  on  the  value  of  u^  To  insure  that  this  initial  condition  is  satisfied 
by  the  approximation,  the  error,  A(uv),  between  the  Fourier  solution  and  the  initial 

condition  at  x  =  0  is  calculated  and  added  to  the  total  sums  for  all  x  in  [0,1]. 

Examples  of  the  results  obtained  for  several  values  of  uv  are  shown  and  are 

compared  with  numerical  results  for  the  same  values  in  figure  3.9.  This  method  yields 
reasonably  good  approximations  to  the  zero-order  solution  and  retains  the  main  features 
which  are  dependent  on  uv.  But  when  this  method  was  used  to  represent  Uq  in  the  first- 

order  analysis,  the  accumulated  errors  were  excessive.  Thus  to  retain  accuracy,  only 
numerical  methods  were  used  in  the  following  first-order  analysis. 

3.5  Time  Scale  Dependence 

Throughout  this  analysis  the  zero-order  equation  has  been  assumed  to  be 
independent  of  the  laboratory  time,  t,.  The  object  of  this  section  is  to  investigate  this 

assumption  by  using  the  coupled  time  scale  terms  of  the  first-order  equation.  If  there  is  a 
possibility  that  there  may  be  a  dependence  upon  ts  we  must  uncover  it.  Beginning  from  the 
elliptical  small  oscillation  result  for  Uq,  a  generalized  2-scale  analysis  is  applied  such  that 

uo  =  uo^.x;  Uv). 

3.5.1  Small  Oscillations 

One  begins  the  argument  by  returning  to  equation  (3.22)  and  differentiating  the 
whole  equation  partially  with  respect  to  x.  This  leaves  the  right  side  of  the  equality  equal  to 
zero.  Integrating  this  again  leaves  a  constant  of  integration  which  is  at  most  a  function  of 
t,,  such  that 


Normalized  bubble  radius.  uQ(x). 
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Numerical  and  Fourier  series  approximations  of  uQ. 


0.0  0.2  0.4  0.6  0.8  1.0 

Normalized  bubble  time,  x. 


Figure  3.9  Comparison  of  Numerical  Zero-Order  Solutions  with  50-Term 
Truncated  Fourier  Series  Approximations  for  Several  uv  Values. 


58 


_  +  j^°  —  CJ  =  1  +  A(ts)  , 


(3.46) 


where  A(o)  =  0  to  hold  with  the  initial  conditions.  Separating  this  into  two  integrals  as  in 
the  previous  analysis  and  integrating,  one  finds  the  bubble  time  plus  another  constant  of 
integration  depending  on  the  laboratory  time,  such  that 


f  *  =  sin'1  £  +  B(Q  +  | 
where  B(0)  =  0  from  the  initial  conditions  and  now 


(3.47) 


*  =  ”o  ~  c- . 

aVl  +  A(ts) 

Substituting  for  and  solving  for  uq,  one  finds  the  2-scale  representation  as 


(3.48) 


u0(ts,t;  uv)  —  c  -  a^l  +  A(ts) 


\™(r 


B(ts) 


(3.49) 


In  order  to  suppress  secular  terms  in  the  zero-order  solution  the  coupled  time  scale  terms 
in  the  first-order  expansion  equation  must  be  set  equal  to  zero.  To  ensure  this,  they  are  set 
equal  to  zero  here  and  evaluated  using  the  generalized  2-scale  equation  (3.49)  for  u0. 

The  coupled  time  scale  terms  from  the  first-order  expansion  equation  set  to  zero  are 


2u»S  +  3§^  =  0' 


(3.50) 


Substituting  equation  (3.49)  and  its  partial  derivatives  into  equation  (3.50),  we  find  that 


be  TdA  .  f'b 
TTaIA 


x  -  B  1-2(1  +  A)^cosl^x  -  B 


5  .  dA  .fb„  tj 

“2abdfsu,lit  '  B 


x  -  B 


ab(l  +  A) 


fM- 


2  cost  -  x  -  B 


j  =  °. 


(3.51) 
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Expanding  the  arguments  using  trigonometric  identities  and  factoring  out  terms  involving 
the  bubble  time,  x,  we  can  write  the  above  equation  as 

be  ffdA  „  ,  dB  .  .  b 

V— •  aIIa  cosB  '  2(1  +  A>d£s,nB|smiT 

-  [a$ sin B  ■  2(1  +  A)  3^  cosBj  cosa  TJ 
•  'i  ab{[a^cos2B]  •  2<>  +  A)S|sin2B  sin2jx) 

-  [s?  sin2B  '  2(1  +  A)  ^  cos2b1cos2^ x 

-[t(1  +  A>a?]  -  °-  (352) 

The  trivial  solution  to  this  equation  is  dA/dtj  =  dB/dts  =  0.  If  a  more  general 
solution  exists,  it  must  exist  for  all  time  x  and  ts.  Since  the  constants  in  front  of  the  major 
brackets  are  not  equal  to  zero  and  the  periodic  terms  involving  the  bubble  time  are  linearly 
independent  functions  and  not  all  are  equal  to  zero  at  the  same  time,  the  only  way  to  satisfy 
the  equality  is  to  equate  each  square  bracketed  term  to  zero.  This  leads  to  five  independent 
equations  all  of  which  are  equal  to  zero. 

The  first  two  square  bracketed  expressions  lead  to  two  equations  which  can  be 
solved  for  the  two  unknowns,  A  and  B.  Looking  at  the  first  of  the  square  bracketed 
expressions,  we  see  that  it  can  be  written  as 

1  dA  „dB  sinB  _ 

(1  +  A)dts  '  ^dt,  cosB  "  u  ’  (3.53) 

An  integration  of  (3.53)  with  respect  to  laboratory  time  yields 

In  (1  +  A)  +  In  (cos2  B)  =  C  ,  (3.54) 

where  C  is  an  integrating  constant  Evaluating  this  for  the  initial  values  of  A  and  B,  one 
sees  that  C  =  0.  Solving  this  equation  for  A  in  terms  of  B,  one  has 
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A  = 


sin2B 


cos2B 

Its  derivative  with  respect  to  laboratory  time  is  then 

dA  _  dB  sinB 


(3.55) 


dtj  ~  ‘‘'dtj  cos3B  (3-56) 

These  two  expressions  can  now  be  substituted  into  the  second  square  bracketed  expression 
from  equation  (3.53)  to  obtain  an  equation  solely  in  terms  of  B.  It  is 


_dB  sin2B  sinB  _  _ 
dts  cos3B  +  cos2B. 


(3.57) 


This  equation  is  satisfied  if  either  dB/dtj  =  0,  or  the  part  in  brackets  is  equal  to  zero. 
Equating  the  bracketed  term  to  zero  we  get 


tanB  =  -  1  , 


which  leads  to 


B  =  (3  +  4n)  j  and -(l  +  4n)*,  n  =  0,  1,  2,  3, ... 

But  these  values  of  B  violate  the  initial  condition  which  requires  that  B(0) 
the  only  solution  to  equation  (3.57)  for  all  tg  is 

dB=0 

dts  "  0  ’ 

or,  from  the  initial  value, 


(3.58) 

(3.59) 

0.  Therefore, 

(3.60) 


B(g  =  0  .  (3.61) 

Substituting  this  back  into  equation  (3.55)  one  gets  the  constant  value  of  A  for  all  ts  to  be 


A(t,)  =  0  ,  (3.62) 

which  satisfies  the  initial  value. 

It  now  must  be  determined  whether  the  other  three  expressions  derived  from 
equation  (3.52)  hold  for  the  same  values  of  A  and  B.  Applying  the  previous  procedure  on 
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the  third  and  fourth  square  bracketed  expressions  from  (3 .52),  first  find  the  expression  for 
A  to  be 

,  1  -  cos2B 

A  “  cos2B  •  (3.63) 

Substituting  this  and  its  derivative  into  the  fourth  expression  we  get  for  B 

2®s^2B=0.  (3.64) 

This  equation  says  that  since  secant  squared  never  equals  zero  for  any  B,  again 
dB/dts  =  0,  and 


b  eg  =  0  . 


Putting  this  back  into  equation  (3.63)  one  again  gets  the  result 


(3.65) 


A(ts)  =  0  . 


(3.66) 


It  is  also  recognized  that  the  fifth  bracketed  expression  from  equation  (3.52)  is  also  satisfied 
for  both  these  conditions.  Applying  these  results  to  equation  (3.49)  one  gets  the  same 
expression  derived  earlier  which  was  independent  of  the  laboratoiy  time,  V  Evidently  in 

the  small  oscillation  case,  the  zero-order  solution  is  indeed  independent  of  the  laboratory 
time,  4.  It  now  remains  is  to  determine  whether  or  not  this  is  the  case  for  the  formal  zero- 

order  equation  over  the  full  range  of  the  parameter,  uv. 

3.5.2  The  Integrabilitv  Condition:  Basic  Structure  and  Most  General  Solution 

The  integrability  condition  is  usually  satisfied  by  substitution  of  an  explicit  zero 
order  solution,  Uo(t,ts)  in  order  to  find  the  dependence  of  Uq  upon  tj.  In  the  present  case 
this  can  be  done  only  for  values  of  u*  very  near  1  where  the  small  oscillation  solution 
applies.  As  preparation  for  the  determination  the  form  of  Uq  with  respect  to  ts  for  all 
values  of  Uy  we  now  consider  the  properties  of  Uq  that  follow  from  the  partial  differential 
equation  itself,  namely: 


2u o 


yuo 

dxdt. 


+  3 


duo<too_  n 

*  dt,"U- 


(3.50) 
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If  Uq  should  happen  to  be  independent  of  t,  this  integrability  condition  is  certainly 
satisfied.  But  are  there  also  forms  of  Uo(t,tg)  depending  on  tj  as  well  as  x  which 
satisfy  the  integrability  condition  ?  That  is  the  question  to  be  answered  and  some 
knowledge  about  the  basic  properties  of  the  partial  differential  equation  and  the  forms 
of  the  integrability  functions  that  will  satisfy  it  may  be  helpful. 

This  partial  differential  equation  is  a  hyperbolic  quasilinear  equation  in  normal 

form  with  characteristic  directions  on  the  lines  x=  constant  and  ^  =  constant.  In  order 

9u0  3uo  d^UQ  92uo  92uo 

to  see  this  we  let  p=  .  q  =  ^  ,r=-^  ,s=  and  w  =  . 

Then  the  two  equations, 


and 


dp  =  r  dx  +  s  dtg 


(3.67) 


dq  =  sdt  +  wdt„  (3.68) 

result  directly  from  the  chain  rule  applied  to  the  formulas  for  p  and  q.  The  partial 
differential  equation  can  now  be  expressed  as 


s  +  ^P0  =  0' 

Suppose  next  that  one  has  a  curve  T  in  the  t,t,  plane 

duo 

along  which  uq  and  its  normal  derivative,  ,  are 


prescribed.  This  is  the  same  as  saying  that  uq  ,  p  and  q  as 
well  as  dp  and  dq  are  known  along  T.  Consequently,  along 
T  we  have  three  equations  in  three  unknowns: 


(3.69) 


0  +  r  dt  +  s  dt  *  dp  , 

(3.70) 

0  +  s  dt  +  wdt  =  dq  , 

(3.71) 

o  +  S  +o 

(3.72) 
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These  equations  enable  us  to  determine  the  derivatives  r,  s  and  w  everywhere  along  T, 
provided  that  A,  the  determinant  of  these  equations,  never  vanishes  along  T.  When  this 
is  true  T  is  called  an  ordinary  curve  and  it  is  possible  to  find  derivatives  of  u0  of  all 
orders  along  T  and  one  can  therefore  expand  Uq  as  a  Taylor  series  about  every  point  on 
T  in  order  to  determine  the  solution  on  a  strip  in  the  x.t,  plane.  The  width  of  this  strip  is 
determined  by  the  radius  of  convergence  of  the  series  at  every  point  along  T.  Character¬ 
istic  directions  at  any  point  on  T  are  defined  as  those  directions  for  which  A  vanishes.  In 
such  a  situation  the  derivatives  of  Uq  can  not  be  found  and  the  computation  of  the  solu¬ 
tion  as  described  above  can  not  be  carried  out.  In  the  present  case  if  any  point  on  T  is 
tangent  to  such  a  direction  we  must  have 

dt  dts  0 

A  =  0  dx  dts  =  dx dt,  =  0,  (3.73) 

0  1  0 


and  the  solution  can  not  be  continued  as  indicated  above.  Clearly,  all  lines  x  =  constant 
and  t,  *  constant  must  be  characteristic  lines  as  stated  at  the  outset. 

When  r  is  an  ordinary  curve  one  can  use  Cramer's  rule  in  order  to  find  the 
value*,  of  r,  s,  and  w  at  any  point  on  T.  Thus  we  have  for  the  present  problem  the 


solution. 


(A)  s  = 


dx  dp  0 
0  dq  d  t. 


(3.74) 


On  a  characteristic,  since  A  =  0,  the  best  that  we  can  hope  for  is  that  solution  for  s  will 
be  simply  indeterminate  on  T.  Such  will  be  the  case  if  we  also  insist  that 
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dx  dp  0 

0  dq  d  tj 

0  -~~pq  0 

2uo 

along  the  characteristics.  This  characteristic  relationship  provides  the  key  for  the  deter¬ 
mination  of  the  dependence  of  Uo  upon  the  laboratory  time  tg.  In  order  to  see  how,  we 
recall  that  all  derivatives  in  the  zero-order  partial  differential  equation  determining  u0  are 
taken  with  respect  to  t.  Therefore  solution  Uq  is  reduced  to  a  quadrature  by  means  of 
two  partial  integrations  with  respect  to  x  and  its  only  dependence  on  ts  is  found  in  the 
constants  of  integration  A(ts)  and  B(ts).  But  in  terms  of  the  integrability  condition  from 
the  first  order  perturbation  equation,  a  partial  integration  over  x  implies  that  d^  =  0.  But 

d  u0  9uo 

d  u0  =  -JY  d  x  +  j^~dts=  pdx  +  qdtj. 

Therefore  we  have 


qdts  =  duo  -  pdx  . 
Therefore  the  characteristic  relationship, 

dTdt*  2^Pq  =0’ 

becomes 


dx  — p(duo  -  pdx)  =0. 

3 

But  the  factor  dx  - —  p  is  certainly  not  generally  zero.  Therefore  the  condition 

2Uo 

(duo  -  pdx)  =0  suggests  that  Uq  is  independent  of  t,  in  general. 


The  fact  that  the  integrability  condition  of  Equation  (3.50)  has  a  rather  symmetric 

form  and  has  a  cross  second  derivative  for  its  highest  order  term  suggests  that  one  might 

find  a  general  form  for  u^x.t,)  which  is  analogous  to  d'Alembert's  well  known  solution 

duo 

for  the  wave  equation.  In  order  to  exploit  this  symmetry  we  can  let  v  =  as  before. 
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(3.76) 


Then  the  integrability  condition  is 

V  9ts  Uo  dts 

0 

Evidently  this  equation  can  be  rewritten  as  r-  In  (v  u3/2)  =  0,  from  which  one  finds  that 

dtj 

=  f(x).  A  second  partial  integration  with  respect  to  x  leads  to 


u05/2  =  J  2  f(x)  dx  +  =  g(x)  +  h(g. 

Consequendy,  the  integrability  condition  wili  be  satisfied  by  any  functions  g(x)  and  h(ts) 
in  an  equation  of  the  general  form, 


uo(t,g  =  [g(x)  +  h(g]2/5.  (3.77) 

That  this  result  is  the  most  general  form  of  the  integrability  function  can  be  verified 
direcdy  by  differentiation  to  show  that 

2u°  = -^  [gW  +  h(g] -6/5  and  that  ^  [gW  +  h(g] 

As  we  have  noted,  equation  (3.77)  is  expressed  in  characteristic  coordinates.  In 
order  to  exhibit  a  form  which  gives  the  solution  on  an  ordinary  curve,  Rxj.ti)  as  illus¬ 
trated  above,  we  need  a  transformation  between  the  characteristic  and  T  coordinates, 
x  =  Mxj.tj)  and  tj  =  yfx^t!).  For  example,  suppose  we  take  x  =  (Xj+  tx)/V  2  and 
h  ~  (xf  liW  2 ,  corresponding  to  T  being  a  line  at  45°  with  respect  to  the  horizontal 
and  upon  which  q  =  0.  Then  one  can  write 


UoCti.ti)  =  [gidi+  h)  +  h^Xj-  ti)]2/5, 

where  the  suffixes  on  gj  and  hj  denote  that  the  T  coordinates  have  been  rescaled  in  order 
to  cause  uo(Xi,tj)  to  resemble  more  closely  the  classical  d'Alembert  solution.  This  formal 
resemblance  must  not  be  allowed  to  obscure  the  fact  that  physically,  the  proper  coor¬ 
dinates  are  x  and  t,  and  that  the  question  of  whether  or  not  the  zero-order  solution 
depends  on  t,  as  well  as  x  is  uppermost  in  our  minds. 
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Two  successive  partial  integrations  of  Equation  (3.1)  with  respect  to  x  have  given 
the  quadrature, 

JSP 

x  +  b(ts)  =  (1  +  q)3/4  — .  -S*2  d5  (3.78) 

'V  2uJln^-^+i  +  a(ts) 

J  ' 

l 

which  is  certainly  not  of  the  general  form  required  by  the  integrability  function.  Equation 
(3.78)  can  be  brought  into  an  acceptable  form  if  a(tj)  =  b(tj)  =  0.  This  choice  satisfies 
the  initial  conditions  and  we  may  then  write  the  solution  symbolically  as  u0(ts,x)  =  G(x), 
where  it  is  supposed  that  the  integral  can  be  evaluated  and  inverted  (numerically  if  need 
be.)  More  generally  one  could  suppose,  because  of  the  small  oscillation  result  discussed 
above,  that  this  inversion  might  have  a  form  such  as, 

u0(x,ts)  =  A(g  F[x  +  B(tg)]. 

The  general  form  of  Equation  (3.77),  appears  to  preclude  this  possibility,  however;  and 
this  conclusion  is  reinforced  by  the  direct  use  of  this  tentative  form  of  u0  in  the 
integrability  condition,  equation  (3.50). 

3.5.3  A  Direct  Approach  to  the  General  Solution 

Since  the  inverse  form,  x[uv;  Uq,  a(ts),b(ts)],  has  been  reduced  to  the  quadrature 
of  equation  (3.78),  one  can  use  the  Leibnitz  rule  in  order  to  test  further  the  finding  ob¬ 
tained  thus  far  in  the  hope  that  if  there  is  a  dependence  on  tg ,  it  might  at  least  be  found 
numerically.  We  shall  start  our  direct  approach  using  the  closed-form  result,  v  =  9uq/9 x, 
which  we  shall  now  write  as 

(1+  q)3/2  v2  =  [2Uylnuo  -  uj  +  1  +  a(tg)]/u30  . 


Therefore  the  mixed  second  derivative  is  found  to  be 


(3.79) 
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2  2 
Uv“U0 


3v  >9uq 


1  da 


vuq  (1+  q)3/2  (1+  q)3/2  2uqV  dts 


so  that  the  second-order  term  in  the  integrability  condition  is 

2(uJ-uJ) 

92u0  _ _  9uo  1  1  da 

vuja+q)**  *.  +a+q)3/2„Jvd>s' 


1  da 

(3-80) 


The  second  term  in  the  integrability  condition,  involving  the  product  of  duo/dr  and 

9ua  8u0 

9uo/9ts,  can  be  written  as  3  v-^-  which  eliminates  the  term,  -3v^“  ,  in  the  mixed 

second  order  derivative  of  equation  (3.80)  from  the  integrability  condition.  Consequently 

dv  dun 

the  integrability  condition,  when  written  as  2uq^  +  3v^“  =  0,  becomes 


2(uJ-Uq) 


1  1  da 


v  Uq  (1+  q)3/2  dtj  (1+  q)3/2  u2y  dts 


This  result  simplifies  to 


2  2,duo  da  _ 

^“v-V^+oo^  =°- 


(3.81) 


The  only  factor  depending  on  b(tj)  results  from  the  derivative,  dujdts.  After 

applying  Leibnitz's  rule  [12]  to  equation  (3.78)  and  rearranging  the  result,  one  finds  that 

“o 

db  v  da  ff  £  l3/2.* 

^st  =  v_-  +  _(i+  q)3M  —  — - - a -  d^  ,  or  that 

dts  2  dts  auJln^-^+i  +  a^ 


“0 

9up  _  db  v  da  f_d£_ 

at,  -  vdt,  +  2(1+ q)3  dt,J43v3  • 


(3.82) 
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Equation  (3.82)  can  now  be  substituted  into  the  integrability  condition  of  (3.81)  After 
separating  the  terms  involving  the  derivatives  of  a(tg)  and  b(ts),  Then  one  finds  that 

/  'N 


I  (<~uo) 


(l+q)3  J[^v(4,a)P 


>  da  _ ,  2  2.  db  _ 

^+2(nv-Uo)v^  =  0. 


(3.83) 


The  integrability  condition  of  equation  (3.83)  is  a  single  condition  for  the  two  unknown 
functions,  a(ts)  and  b(ts).  If  a  and  b  are  truly  independent  arbitrary  functions  then  they 
must  vanish  identically.  Otherwise,  equation  (3.83)  implies  that  a  and  b  are 
interdependent . 

Suppose  therefore,  that  we  were  to  consider  a  =  a(b).  The  differential  dtj  is  then 
eliminated  from  Equation  (3.83).  If  in  general  the  differentials  da  and  db  are  not  zero 
simultaneously,  equation  (3.83)  will  be  satisfied  only  if  the  coefficients  of  da  and  db 
vanish.  Then  we  would  have  two  equations  for  the  unknowns  a(t,)  and  Uq.  If  these  equa¬ 
tions  were  solved  for  values  of  the  unknowns  near  their  starting  values,  a  =  0  and  u0  =  1. 
Then  one  would  evaluate  Equation  (3.78)  from  1  to  u©  and  so  calculate  the  sum,  tj  +  b(ts), 
with  the  help  of  the  condition  tj  =  ex.  But  he  could  not  determine  %  or  b  without  an  added 
constraint  which  does  not  appear  to  exist  for  this  formulation.  Thus  we  see  that  if  aftg)  and 
b(tj)  are  not  independent  functions  of  integration  ,  the  system  of  equations  is  not  closed. 
Moreover  even  if  the  system  were  closed,  these  two  equations  do  not  vanish  at  the  initial 
condition,  where  a(0)  =  v(0,0)  =  0  and  Uq(0,0)  =  1.  Then  the  requirement  that  the  coeffi¬ 
cients  of  da  and  db  vanish  for  all  t,  and  t  is  certainly  not  satisfied. 

Suppose  next,  that  one  considers  the  differential  equation  (3.83)  and  the  quadrature 

of  Equation  (3.78)  as  two  of  a  governing  system  of  equations.  How  might  the  system  be 

closed?  One  possibility  may  be  suggested  by  the  behavior  of  the  zero- order  radial  velocity, 
duo 

— ,  which  from  the  chain  rule  and  the  condition  dtg  =  edx,  is 
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duo  duo  duo 

dT=  aT  +  Edt^  • 


(3.84) 


The  two  partial  derivatives  of  u<)  in  equation  (3.84)  are  given  by  the  square  root  of 


equation  (3.79)  and  by  equation  (3.82),  respectively.  Equation  (3.84)  can  be  used  to 

determine  the  maximum  and  minimum  values  of  Uq.  These  values  are  given  by  the  two 
duo 

roots  of  =  0,  namely: 


[2u^  lnu0  -  Uq  + 1  +  a(t)]/i^ 
(1+q)3'4 


db 

+  e  vr  + 


dts  2(1 


1 _  da  fd£_  _  (3.85) 

+  q)3  dtg  J  £3  v3  u- 


If  the  roots  of  (3.85)  can  be  found  starting  at  the  initial  point,  the  value  of  umax, 
say,  permits  one  to  assign  a  series  of  equally  spaced  points  over  the  interval  [l,umax]. 
Let  the  first  new  value  of  uq  be  1+h,  where  h  is  the  step  size.  Of  course,  u^*  and  umjn 
must  be  found  at  each  step  because  the  values  of  a  and  b  will  change  with  each  compu¬ 
tational  step.  Therefore  it  is  possible  that  the  size  of  each  step  should  be  found  from 
integration  of  equation  (3.84)  and  equation  (3.85)  should  be  reserved  for  the  determi¬ 
nation  of  Uquh  and  u^.  We  note  in  equation  (3.85)  that  a(t$)  and  b(tg)  are  the 


unknowns.  The  same  is  true  of  equation  (3.83).  Equation  (3.78)  has  a,  b  and  x  as 
unknowns,  provided  that  Uq  can  be  found  from  integrating  equation  (3.83).  Then  by 
making  use  of  ts  =  ex,  we  see  that  the  fundamental  unknowns  are  a,  b  and  x  at  each  step 
and  there  are  three  equations.  Therefore,  the  system  is  closed;  but  it  remains  to  design  an 
algorithm  for  a  step  by  step  solution.  All  of  this  is  conjectural.  Therefore  the  present 
study  is  restricted  to  the  case  in  which  a  =  b  =  0  and  the  integrability  condition  is  indeed 
satisfied.  This  approach  may  be  useful  preparation  for  an  attack  on  this  generalized 
problem,  should  it  be  necessary.  The  present  special  formulation  will  reveal  some  aspects 
of  the  solution  properties  that  must  also  be  dealt  with  even  in  a  generalized  approach. 
This  will  be  especially  true  in  the  neighborhood  of  of  the  initial  point  x  =  tj  =  0  because 


a  and  b  will  be  very  small  in  this  neighborhood. 
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CHAPTER 4 

FIRST-ORDER  EQUATION  ANALYSIS 


The  two-scale  perturbation  expansion  in  section  2.4.2  gives  the  l^rst-order  linear  pertur¬ 
bation  partial  differential  equation  for  u^x.tj)  where  Uj  is  the  first  order  normalized 
dimensionless  bubble  radius,  Uj.  The  two  dimensionless  time  scales  are  ts,  and  x  and  the 
first  order  equation  was  found  to  be. 


d2Uj  3  9u0  fduA  U!  92u0  1  ^uv 

dX2  +  U0  dX  [dx  J+U0  Bx2  +u2(1  +  q)3/2  u2 


=  -2 


32u0  3  j_  9uq  9u0  _ F(t 

dxdts  +  2  u0  dx  +  u0(l  + 


F(ts) 

;i  +  q)3/2  ’ 


with  initial  conditions  of 


Ui(0,0)  =  0 


9u1(0,0)  9ug(0,0)  _  l 

dx  +  c*tl  “  8 


P  =  iwrV(l+q)(l+K)(-  ^]c  =  k  .  (4.3) 


In  equation  (4.1),  F(t*)  is  the  pressure  forcing  function  of  section  2.2  and  the  term 
to  the  right  of  the  equality  in  (4.3)  is  the  flaccid  bubble  response,  from  section  2.3.4,  to 
F(ts) .  The  first-order  equation  depends  on  the  zero-order  equation  and  its  derivatives  as 
well  the  affix  of  the  vortex  point,  uv. 


To  write  this  equation  in  terms  of  the  normalized  bubble  time,  x,  one  applies 
equation  (3.17)  and  its  derivatives  to  equation  (4.1).  In  section  3.5  we  have  taken  the  zero- 
order  solution  to  be  independent  of  the  slow  laboratory  time,  t,,  so  that  those  terms 
containing  derivatives  of  Uq  with  respect  to  t,  must  vanish.  Moreover,  the  mixed  derivative 
terms  in  the  second-order  equation  of  u2  also  vanish  for  die  same  reason.  Thus,  the  first- 
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order  equation  may  be  considered  to  be  independent  of  the  laboratory  time,  ts,  and  can  be 
written  as  an  ordinary  differential  equation  depending  only  on  the  bubble  time,  x.  Applying 
this  and  the  period  parameter  from  equation  (3.20),  we  can  write  a  normalized  form  of  the 
first-order  equation  as 


d2ut  3_  d^o  fduA 
dx2  u0  dx  |^dx  J 


Uj_  d2u0 
Uq  dx2 


+ 


with  corresponding  initial  conditions  of 


1 


X2F(e,x) 


uo 


(4.4) 


Ul(0)  =  0 


(4.5) 


and 


du^O) 

dx 


^Wr(l+q)5/4V(TTK) 


(4.6) 


This  leaves  a  single-degree-of-ffeedom  system  governed  by  a  linear  second-order 
nonhomogeneous  ordinary  differential  equation.  Due  to  the  fact  that  the  zero-order 
equation  has  periodic  solutions,  this  ordinary  differential  equation  has  time  dependent 
periodic  coefficients  that  can  produce  parametric  excitations  in  the  system.  This  type  of 
system  is  discussed  by  Nayfeh  and  Mook  [7]  among  others,  and  as  they  suggest  it  can  be 
solved  using  ihe  Floquet  theoiy  for  the  homogeneous  solution  and  variation  of  parameters 
for  the  particular  solution  [12]. 


4.2  Homogeneous  Solution  from  the  Floquet  Method 
4.2.1  Standard  Form  of  the  Equation 

The  linearity  of  equation  (4.4)  permits  the  forcing  term  to  be  set  to  zero,  leaving  a 
homogeneous  equation  involving  only  derivatives  with  respect  to  the  normalized  bubble 
time,  x.  The  homogeneous  part  of  Equation  (4.4)  can  then  be  written  as 

d2u  .  .  du  ^  ,  .  n 

^2  +  PiW  +  P2(x)u  =  0  ♦  (4.7) 


where 
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The  parameter  p(x)  is  also  a  periodic  function,  dependent  on  the  zero-order  solution 
and  uv.  Figure  4.1  shows  the  form  of  p(x)  for  various  values  of  Uy.  It  is  seen  that  as  the 
value  of  Uy  increases,  portions  of  the  function  p(x)  become  increasingly  steeper  and 

steeper.  This  phenomenon  causes  equation  (4.1 1)  to  become  increasingly  stiff  with  rising 
Uy  values,  and  some  care  must  be  taken  if  this  equation  is  to  be  evaluated  numerically. 

This  standardized  form  of  the  first-order  equation  in  (4. 1 1)  is  among  the  class 
known  as  Hill's  equations,  and  can  have  stable  or  unstable  solutions  depending  on  the 
value  of  Uy  in  p(x).  For  the  special  case  when  p(x)  has  die  form 


Coefficient  p(x)  in  first-order  equation. 


Derivatives,  dp/dx. 
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(4.13) 


p(x)  =  d*  +  2e*cos(2x)  , 

where  d*  and  e*  are  stability  parameters,  equation  (4.11)  reduces  to  a  form  of  Mathieu's 
equation  dealing  with  the  problem  of  vibrations  of  an  elliptic  membrane.  With  the  homoge¬ 
neous  first-order  equation  in  this  standard  form  its  stability  and  its  solution  are  well  known 
and  are  given  by  Abramowitz  and  Stegun  [9]  among  others. 

4.2.2  Small  Oscillation  Stability  Analysis 

As  discussed  in  Nayfeh  and  Mook  [8],  the  stability  of  equations  (4.7)  and  (4.1 1)  is 
the  same.  So  the  stability  of  one  will  be  the  same  as  that  of  the  other.  Using  the  small 
oscillation  results  from  the  zero-order  solution  for  the  case 


uv  =  1  +  8  ,  (4.14) 

where  5  «  1,  a  small  oscillation  form  of  the  standard  equation  can  be  obtained.  Sub¬ 
stituting  the  small  oscillation  expressions  for  Uq  and  1  from  equations  (3.34)  and  (3.37)  into 


equation  (4.1 1)  and  keeping  terms  to  the  second  order  in  8,  one  finds  a  small  oscillation 


form  of  the  homogeneous  first-order  equation  to  be 
^|  +  4jc2£i+8 ycos(2jtx)  +  S2(y-  y|cos(2jtx)  +  ycos(4jtx))J  z  =  0. 


(4.15) 


Equation  (4.15)  can  be  recognized  as  a  Hill's  equation.By  letting  ft  =  tcx  and 
ruling  out  terms  containing  46;  one  can  write  this  equation  to  the  second-order  in  8,  except 
for  the  neglected  46  term,  as 


d2z 

d62 


+  4[l+8*f  +  8  (y  -  d§)cos(2e)]z  -  0  .  (4.i6) 


This  is  a  standard  Mathieu  equation  with  the  form  of  the  coefficients  as  described  in 
equation  (4.13).  The  term  S2  g-cos(4  n  x)  from  equation  (4.15)  was  left  out  of  the 

expression  in  (4.16)  in  order  to  produce  the  proper  form  for  a  Mathieu  equation.  This 
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omission  seems  acceptable  since  similar  results  to  the  present  analysis  were  obtained  while 
only  keeping  terms  up  to  the  first  order  in  8. 

Comparing  the  coefficient  in  equation  (4.16)  to  the  Mathieu  form  in  equation 
(4.13),  we  find  that  the  stability  parameters  are 


8*  =  4(l  +  82f-) 


(4.17) 


e‘=8(n-8f)  . 


(4.18) 


Therefore  both  8*  and  e*  are  functions  of  8. 


The  stability  of  the  Mathieu  equation  can  be  visualized  by  plotting  8*  against  e*.  In 
the  8*,e*  plane,  for  the  present  case  there  are  two  curves  of  neutral  stability  [9]  given  by, 


5*  ■ 4  -  h<e*>2  +  +  - 


(4.19) 


8*  =  4  +  12(e*)2  -  13321  (e*)4  +  •••  * 


(4.20) 


In  the  region  between  these  two  curves  the  solution  to  the  Mathieu  equation  is 
unstable,  while  outside  these  curves  the  solution  is  stable.  Remembering  that 


5  »  uv  -  1  ,  (4.21) 

one  can  plot  equations  (4.17)  and  (4.18)  as  in  figure  4.2  for  various  small  values  of  tv 

Note  that  two  such  curves  are  plotted  in  figure  4.2 .  One  curve  is  shown  as  a  continuous 
curve  with  circled  points  showing  various  uv  values  in  the  interval  [1,1.013].  The  second 
curve  is  shown  as  a  sequence  of  small  black  squares  and  this  curve  covers  the  Uy  interval 
[1, 0.987] .  Included  in  this  figure  are  the  neutral  stability  curves  from  equations  (4.19) 
and  (4.20). 

As8-+  0  and  Uy  -+  1,  the  Mathieu  equation  in  question  yields  neutral  stability 
parameters  of  8*-»  4  as  £*-»  0.  Similarly,  both  neutral  stability  curves  also  approach  a 
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value  of  8*->  4  when  e*-»  0.  This  shows  that  the  homogeneous  part  of  the  first-order 
solution  departs  from  a  point  of  neutral  stability  at  Uy  =  1,  continuing  into  an  unstable 
region  as  uv  increases  or  decreases.  When  Uy  =  1,  the  period  parameter  assumes  the  value 
X  =  tcV2  and  the  Mathieu  equation  in  (4.1 1)  becomes  the  equation  of  a  harmonic  oscillator 
which  is  indeed  neutrally  stable.  It  appears  from  figure  4.2  that  the  small  oscillation 
solution  curves  continue  to  stay  in  the  unstable  region  for  all  higher  values  of  Uy,  although 
the  small  oscillation  approximation  will  break  down  as  the  value  of  Uy  increases  or 
decreases  too  far  outside  of  the  range  [0.98, 1.03]. 

4.2.3  Floauet  Method 

The  homogeneous  first-order  equation  describes  the  behavior  of  a  system  governed 
by  a  linear  ordinary  differential  equation  with  periodic  coefficients.  The  functional  behav¬ 
ior  of  single-degree-of-freedom  systems  can  be  characterized  using  the  Floquet  theory. 

This  theory  determines  general  properties  of  the  solutions  from  systems  governed  by 
equations  of  the  form  (4.7),  without  actually  solving  them.  These  properties  can  then  be 
used  to  determine  numerical  approximations  to  the  solutions.  These  approximations  have 
factors  that  are  periodic,  allowing  them  to  be  used  over  many  cycles,  whereas  other 
numerical  solution  techniques  tend  to  accumulate  excessive  errors  when  used  over  an 
extended  range.  The  nonperiodic  factors  are  exponentials  which  also  permit  accurate 
numerical  evaluation.  The  following  development  is  given  by  Nayfeh  and  Mook  [8],  and 
applied  to  die  equation  at  hand. 

Since  equation  (4.7)  is  a  second-order  linear  homogeneous  equation,  it  has  two 
linearly  independent  solutions  Uj(x)  and  u2(x).  Here  die  subscripts  refer  to  the  two 

independent  solutions  to  the  first-order  homogeneous  equation  and  not  to  the  first  or 
second-order  equations  of  the  two  scale  expansion.  Functions  such  as  U]  and  u2  are 

usually  referred  to  as  a  fundamental  set  of  solutions  because  every  solution  is  a  linear 
combination  of  these  two  solutions,  and  because  they  satisfy  die  fundamental  initial 
conditions  of 


uj(0)  =  1  Uj(0)  =  0 

u2(0)  =  0  u2(0)  =  1 
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(4.22) 


Since  the  zero-order  solution  is  periodic,  equations  (4.8)  and  (4.9)  are  also  periodic 
and 


Pi(x)  =  Pi(x  +  T)  ,  (4.23) 

where  T  is  the  period  and  in  the  normalized  sense  equal  to  unity.  It  follows  that  if  Uj(x) 
and  u2(x)  are  any  two  solutions  of  equation  (4.7),  then  Uj(x  +  T)  and  u2(x  +  T)  are  also 
solutions  of  equation  (4.7).  Also,  if  Uj(x)  and  u2(x)  are  linearly  independent,  then 
Uj(x  +  T)  and  u2(x  +  T)  must  be  linearly  dependent  on  Uj(x)  and  u2(x)  because  a  second- 
order  equation  has  only  two  linearly  independent  solutions.  Therefore,  there  exist  two 
constants  for  each  solution  that  do  not  both  vanish  at  the  same  time  such  that 


{Ui(x  +  T)}  =  [A]{Uj(x)}  , 

where  i  =  1  and  2  and  the  matrix  [A]  contains  the  four  constants 


(4.24) 

(4.25) 


Applying  the  fundamental  initial  conditions  from  (4.22)  and  the  normalized 
period  to  equation  (4.24)  and  its  derivative,  we  find  the  constants  to  be 

*“  =  U,0)  *12  *  “>(1)  (4.26) 

a21  *=  u2(l)  a22  =  u2(l) 

Once  Uj(x)  and  u2(x)  are  known  the  a^  constants  can  be  determined. 

Since  equations  (4.7)  and  (4.1 1)  are  both  a  form  of  Hill's  equation  and  have  the 

same  stability  characteristics,  die  above  process  also  holds  true  for  a  fundamental  set  of 

solutions,  Zj(x)  and  z^x).  In  the  actual  numerical  solution  process  it  is  the  fundamental  set 

of  Zj  solutions  which  are  found,  and  they  are  transformed  back  into  the  U]  form  using 


equation  (4.10).  This  is  done  for  convenience,  since  the  small  oscillation  and  stability 
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results  are  more  easily  obtained  from  equation  (4.1 1),  as  will  be  shown  later.  However, 

Uj  and  u2  will  continue  to  be  used  in  the  development  of  the  solution  to  avoid  confusion. 

It  can  also  be  stated  that  there  exist  fundamental  sets  of  solutions  having  the 
property 

Vj(x  +  T)  =  Xj  vj(x)  ,  i  =  1, 2;  (4.27) 

where  Xj  are  the  eigenvalues  of  the  matrix  [A].  These  eigenvalues  may  be  complex. 

Equation  (4.27)  represents  what  are  called  Floquet  normal  solutions,  and  are  related 
linearly  to  Uj(x)  by  a  2  x  2  constant  nonsingular  matrix,  [P],  such  that 

(Ui(x)}  =  [P]  {Vi(x)J  . 

The  substitution  of  equation  (4.28)  into  equation  (4.24)  yields 

(Vj(x  +  T)}  =  [P]1  [A]  [P]  {Vi(x))  , 
or 

(Vj(x  +  T)}  =  [B]  {Vj(x)}  .  (4.30) 

The  matrices  [A]  and  [B]  are  said  to  be  similar  matrices  because  they  have  the  same 
eigenvalues,  that  is 

l[B]  -  Xj[I]l  =  I  [A]  -  Xj[I]  I  =  0  ,  i  =  1.2.  (4.31) 

The  eigenvalues  of  matrix  [A],  which  will  also  satisfy  matrix  [B],  must  now  be  found1 .  It 
follows  from  linear  algebra  that  a  nonsingular  matrix  [P]  can  be  chosen  so  that  [B]  will 
have  its  simplest  possible  or  Jordan  canonical  form.  This  form  depends  on  the  eigenvalues 
of  [A]  which  are  found  by  solving 

1(A)  -Xi(I]l  =  0  ,  i  =  1, 2.  (4.32) 

1  It  is  hoped  there  will  be  no  confusion  between  X  for  the  period  parameter  and  Xj  for  die 
eigenvalues. 


(4.28) 

(4.29) 


A(x)  =  ZjCxJ^x)  -  z1(x)z2(x)  =  constant  .  (4.39) 

Convening  this  back  to  the  Uj  solutions  using  equation  (4.10),  one  finds  that  the  expression 
is  still  equal  to  a  constant  Evaluating  this  and  using  the  initial  conditions  of  (4.22)  one 
gets 


A(x)  *  1  , 


(4.40) 
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which  is  as  it  should  be  in  the  case  of  a  Hill's  equation  [8].  This  is  an  important  result.  It 
can  be  used  later  to  indicate  roughly  the  error  in  the  numerical  solutions  for  Uj  and  u2. 

The  stability  parameter,  a,  is  used  to  determine  the  value  of  the  eigenvalues,  Xj,  of 


[A].  With  A  =  1,  equation  (4.33)  yields 


X12  =  a  ±  V a2  -  1  , 


(4.41) 


where  the  roots  are  related  by 


Xj  X2  =  1  . 

If  a  =  ±  1 1 1,  equation  (4.41)  yields  only  one  eigenvalue,  namely 


(4.42) 


X  =  a  =  ±  1  . 


(4.43) 


If  a  *  1,  equation  (4.41)  yields  two  distinct  eigenvalues,  and  [B]  takes  on  the  diagonal 


=  [o  x°J  • 


(4.44) 


It  follows  from  equation  (4.30)  that 

vj(x  +  nT)  =  X,"  vj(x)  , 

where  n  is  an  integer  and  i  =  1  and  2.  Then  as  x  -»<*>,  n  and 


0  if  |XJ  <  1 

oo  if  |Xj|>  1  * 


(4.45) 


(4.46) 


WhenXi  *  1,  Vj  is  periodic  with  period  T.  When  Xj  =  -1,  v;  is  periodic  with  period 
2T.  Thus  the  cases  Xi  «  X$  =  ±1  separate  stable  from  unstable  solutions  for  the 
fundamental  set  and  are  called  transition  values. 

It  remains  to  express  equation  (4.27)  in  normal  Floquet  form.  First  we  multiply 
the  respective  equations  by  exp[-  y(x  +  T)],  which  yields 


e*  Yi(x  +  T)  y.(x  +  t)  ,  c-  T,x  V.(X)  § 


(4.47) 
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where 


c'YiT  =  X|  (4.48) 

and 

Yi  =  Jp  In  Xj  .  (4.49) 

The  term  y,  in  equation  (4.49)  is  called  the  characteristic  exponent.  The  right  side  of 
equation  (4.47)  is  periodic  with  period  T,  and  can  be  expressed  as 

e'  *»x  Vj(x)  =  4>j(x)  ,  (4.50) 

where  <J>;(x  +  T)  =  <j>j(x).  Hence,  Vj(x)  can  be  expressed  in  the  normal  form 

Vj(x)  =  eY»x  <t>i(x)  .  (4.51) 

For  the  case  when  Xj  =  ±1,  similar  reasoning  yields  the  same  result  for  Vj(x),  but 
v2(x)  will  have  the  form 

v2(x)  =  eY2x^2(x)  +  3£T  ♦»(*)]  (4.52) 

due  to  the  double  idol  These  give  expressions  for  the  approximate  normal  Floquet 
solutions  of  equation  (4.7)  derived  from  the  properties  of  the  numerical  solutions,  Uj(x) 
and  U2(x).  They  are  made  up  of  well  defined  exponential  functions  multiplied  by  periodic 
functions,  fy. 

The  stability  of  the  fundamental  set  of  solutions  in  Floquet  form  from  equations 
(4.51)  and  (4.52)  can  be  analyzed  by  calculating  the  value  of  a  Whenlal  >  l,the 
absolute  value  of  one  of  the  Xj  is  larger  than  one  while  that  of  the  '"her  is  less  than  one, 
according  to  equation  (4.42).  This  leads  to  and  y2  having  equal  magnitudes  but 

opposite  signs.  Thus,  one  of  the  normal  solutions  is  unbounded  and  the  other  is  bounded 
as  x  ->  ®°,  as  evidenced  by  equation  (4.46).  When  la!  <  1,  the  roots  from  equation 

(4.41)  are  complex  conjugates  having  unit  moduli.  As  stated  before,  the  transition  from 
stability  to  instability  occurs  for  lal  ■  1,  which  corresponds  to  the  repeated  roots 
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X,j  -  X.2  =  1  and  A,!  =  =  -1  so  that  Yj  =  Y2  =  Oorirc.  The  former  case  has  a 

periodic  normal  solution  with  period  T,  while  the  latter  case  has  a  periodic  normal  solution 
of  period  2T.  The  locus  of  transition  values  where  lal  =  1  is  a  very  important  outcome  of 

this  method  and  can  be  used  to  determine  critical  flow  values  that  will  produce  neutrally 
stable  solutions.  The  stability  parameter,  a,  the  eigenvalues,  Xj,  and  the  characteristic 
exponents,  y„  can  be  determined  from  equation  (4.1 1)  by  numerically  calculating  the  two 
linearly  independent  solutions,  zl  and  Z2,  during  the  first  period  of  oscillation.  With  the 

help  of  these  solutions  and  their  first  derivatives  at  x  =  T,  one  can  calculate  a  and  A  from 
equations  (4.34)  and  (4.35).  Then  the  X's  and  Ys  can  be  determined  from  equations 
(4.41)  and  (4.49)  respectively. 

The  solutions  ux  and  u2  can  be  found  by  taking  the  numerical  results  for  zx  and  Z2 
from  equation  (4.1 1)  and  transforming  them  into  Uj  and  u2  with  the  help  of  equation 
(4.10).  But,  due  to  numerical  instabilities  in  the  adaptive  Runge-Kutta  algorithm  used 
here,  errors  will  accumulate  in  the  numerical  process  because  the  equation  is  stiff  in  cetain 
limited  x  intervals  duning  the  integration  over  a  period,  [0,1].  The  degree  of  stiffness 
increases  as  uv  departs  more  and  more  from  unity  .This  aspect  of  the  work  will  be  dis¬ 
cussed  further  below.  Over  the  full  range  of  the  forcing  function  the  solutions  need  to  be 
extended  over  many  periods  in  x,  which  will  greatly  reduce  accuracy  as  the  errors  keep 
accumulating.  To  overcome  this  error  accumulation  die  numerical  solutions  can  be 
transformed  into  normal  or  Roquet  solutions  which  in  general  have  periodic  factors  as 
noted  above.  This  would  require  numerical  solutions  of  ut  and  u2  for  only  one  period  of 

oscillation,  thus  minimizing  the  errors.  This  is  also  all  that  is  required  in  the  calculation  of 
the  previous  stability  parameters  and  characteristic  exponents. 

It  then  remains  to  find  the  components  of  the  Roquet  solutions  in  equation  (4.51) 
in  terms  of  the  numerical  Ui  solutions.  The  Ys  have  already  been  developed.  Expressions 
for  die  periodic  functions,^,  must  then  be  found.  This  is  begun  by  relating  the  uj 
solutions  to  the  V;  solutions  using  equation  (4.28).  To  find  the  components  of  the  matrix 
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[P]  in  (4.28),  the  eigenvectors  corresponding  to  the  eigenvalues  of  matrix  [A]  are  found, 
such  that 


IX[IJ  -  [A] I  [P]  =  0  , 


(4.53) 


where,  for  cases  i  =  1  and  2 


[P] 


-  fPlil 

LP2lI  ' 


(4.54) 


Expanding  equation  (4.53)  for  both  values  of  i,  four  equations  are  obtained  and  can  be 
expressed  as 

a  1  1 P 1  i  +  al  2P 2i  =  ^iPli  [ 

a21  P  li  +  a22P2i  =  ^iP2i  j 
From  this  set  of  equations,  we  find  the  ratio 


(4.55) 


Pli 


a12 


P2i  ‘  al  1 
This  leads  to  the  modal  matrix 

[P]  = 


'  a22 
a2I 


(4.56) 


a12  a12 

Xi  -  a22  ^2  ‘  a2 2 

-  a1  j  X2  "  ai  1  . 

a21  a21 

(4.57) 


Either  of  these  matrices  can  be  used  to  transform  between  variables,  or  their  average  could 
be  used  to  tty  and  minimize  errors  in  a  computer  code.  Only  the  first  form  is  used  here. 

From  equation  (4.28)  it  is  seen  that  the  inverse  of  [P]  is  needed  to  express  the 
solutions  in  Floquet  form.  The  inverse  of  [P]  is 


,1 


1 

X2  -  ai  1 

”  a12 

1 

a21 

-(X,2 

-  a22) 

AP 

.  '(^-t  *  an) 

a12  . 

“  AP 

.  "  a21 

Xi  - 

a2  2  . 

(4.58) 


where 


AP  =  a12(^2  •  ^-t) =  a2i  (^i  ‘  ^2)  •  (4.59) 


For  convenience,  the  first  of  these  is  chosen  for  use  in  the  following  development 
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Applying  equation  (4.58)  to  equation  (4.28),  one  finds  that  the  Floquet  normal 


solutions  in  terms  of  the  numerical  solutions  are 


vj(x)  = 


v2(x)  = 


(X2  -  an)u!(x)  -  a12u2(x) 
a12 (^2  "  ^l) 

-(Xj  -  an) ut(x)  +  a12u2(x) 

a12(^2  '  ^l) 


Equation  (4.51)  can  be  solved  for  4>;.  It  is  found  to  be 


(4.60) 


(4.61) 


<t>i(x)  =  e'yix  vj(x)  , 


and  its  first  derivative  is 


-  TfiVjCx)]  . 


(4.62) 


(4.63) 


The  normal  Floquet  solutions  and  their  derivatives  in  terms  of  the  numerical  solutions  from 

equations  (4.60)  and  (4.61)  can  then  be  substituted  into  equations  (4.62)  and  (4.63).  Both 
d<b; 

0i  and  are  periodic  functions  in  x  and  we  only  need  u4  numerical  solutions  over  the 

first  period  to  be  continued  indefinitely. 

In  summary,  two  linearly  independent  solutions  to  equation  (4.7)  can  be  found 
numerically,  by  starting  with  initial  conditions  for  a  fundamental  set  and  integrating  over 
one  period  of  oscillation.  The  value  of  these  solutions  when  x  =  T  =  1  can  be  used  to  find 
the  characteristic  parameters  of  the  solutions,  namely:  the  stability  parameter.a ,  the 
eigenvalues,  Xj,  and  the  characteristic  exponents,  Yj.  These  solutions  can  then  be  used  to 
determine  a  periodic  function  which,  along  with  the  characteristic  exponent,  can  be  used  to 
express  an  approximate  normal  Floquet  form  of  the  solutions  which  can  be  continued  over 
many  periods.  It  now  remains  to  find  the  particular  solution  to  the  equation,  to  combine  it 
with  complementary  function  and  to  evaluate  the  constants  of  integration,  A  and  B,  in  the 
complementary  function  in  terms  of  the  initial  conditions. 
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To  find  the  particular  solution  to  the  nonhomogeneous  form  of  the  first-order 
equation,  one  can  use  variation  or  parameters  [12].  Ordinarily  it  is  impossible  to  find  a 
particular  integral  of  a  nonhomogeneous  lint  or  differential  equation  by  inspection.  But 
since  a  form  of  the  homogeneous  solution  can  be  found  by  the  Floquet  theory,  the  general 
variation  of  parameters  method  enables  one  to  replace  the  constants  of  integration  in  the 
complementary  function  with  the  functions  a(x)  and  b(x).  This  is  done  by  assuming  a 
particular  solution  to  the  first-order  equation  of  the  form 


Up(x)  =  a(x)uj(x)  +  b(x)u2(x)  ,  (4.64) 

where  again  uj(x)  and  u2(x)  are  the  two  fundamental  solutions  to  the  homogeneous 


solution,  and  a(x)  and  b(x)  are  functions  of  x  that  satisfy  the  additional  condition 

da  .  db  A 
£“■  +  su2  =  0  . 


(4.65) 


Differentiating  equation  (4.64)  with  respect  to  x  two  times  and  applying  the  condition  in 


equation  (4.65)  one  gets 


duD  dui  ,  du2 

-j  =  a-r1  +  b 
dx  dx  dx 


(4.66) 


d2up  _  d2ut  da  duj  d2u2  db  du2 
dx2  a  dx2  *  dx  dx  dx2  +  dx  dx 

The  particular  solution  must  satisfy  the  original  forced  first-order  equation. 

Substituting  these  results  inm  equation  (4.4)  one  has 


(4.67) 


d2Uj 
1  dx2 


dui  d2u2  du2 

+  Pi(x)  d5f  +  P2(x)u!  +  b  +  Pl(x)-5r  +  p2(x)u2 


daduj  dbduz 

+  dx  dx  +  dx  dx  "u0  F(e,x)  ' 


(4.68) 
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But  both  Uj  and  u2  arc  solutions  of  the  homogeneous  equation,  so  that  the  first  two 
terms  in  equation  (4.68)  are  zero  according  to  equation  (4.7).  Therefore  equation  (4.68) 
is  reduced  to 


da  duj  db  du2 
dx  dx  + 


dx  dx  *  57F<£-x>  • 


(4.69) 


Equations  (4.65)  and  (4.69)  yield  two  equations  for  the  two  unknowns,  da/dx  and 
db/dx.  They  can  be  written  as 


‘  “1  u2  ' 

- 1 

_ 1 

i 

» 

O 

_ l 

ii 

X2 

.  «l  u2  . 

1 - 

—F(e,x) 
Luo  J 

(4-70) 


The  determinant  of  this  system  is 


A  =  Uju2  -  uju2  .  (4.71) 

This  is  the  Wronskian  of  the  homogeneous  equation  from  (4.39)  and  is  equal  to 
unity.  Equation  (4.70)  yields  two  equations.  Integrating  them,  one  finds  a(x)  and  b(x)  to 
be 


0 


and 


b(x)  =  X2 


’  «i(E) 
Uo(5) 

6 


F(e,5)  d5  , 


(4.72) 


where  a(0)  =  b(0)  =  0.  Substituting  these  into  equation  (4.64)  one  gets  the  particular 
solution  of 


uD(x)  =  X4 


-uj(x) 


^  4(5) F(e’^  +  U2(x)  J  4(5)  F(e’^ 


6 


(4.73) 
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When  x  =  0,  both  integrals  in  (4.73)  vanish  leaving 


up(0)  =  0  . 

Differentiating  equation  (4.73)  to  obtain  the  particular  growth  rate,  we  find 


dx 


-Ui(x) 


f  u2(S) 


i 


«o(4) 


F(e£)d5  +  i2(x) 


f  ui(4) 

uq(4) 


F(e£)  d£ 


(4.74) 


(4.75) 


Both  of  these  integrals  also  vanish  when  x  =  0,  thus 

=  0  . 


dur 

dx 


x  =  0 


(4.76) 


Equations  (4.73)  and  (4.75)  represent  the  particular  integral  and  its  derivative  for  the 
normalized  first-order  equation.  They  can  be  combined  with  the  complementary  function 
to  obtain  a  complete  integral  of  the  differential  equation  and  then  the  initial  conditions  can 
be  used  to  solve  for  the  constants  of  integration,  A  and  B. 

To  be  able  to  continue  the  particular  solution  numerically  over  many  cycles  we 
avail  ourselves  of  the  above  noxmal  or  Floquet  solution.  Using  it  we  can  transform  the 
particular  solution  into  normal  or  Floquet  form.  To  do  this,  we  substitute  Uj ,  u2  and  their 

derivatives  into  the  particular  solution.  The  matrix  [P]  from  equation  (4.57)  is  used  to  get 
them  in  terms  of  the  Floquet  variables  v}  and  v2,.  According  to  the  relationship  in 

equation  (4.28), 


uj(x)  =  a12[vj(x)  +  v2(x)] 


(4.77) 


U2(x)  =  (Xj  -  an)v!(x)  +  (fa-  an)v2(x)  .  (4.78) 


Substituting  these  expressions  back  into  equations  (4.73)  and  (4.75),  we  find 
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Up(x)  =  X2a12(Xj  -  X2) 


X 

X 

V^x) 

[JgFfcOdt 

6 

6 

.  (4.79) 


and  its  derivative 


^  =  X2a12(Xj  -X2) 

1 

<• 

\ 

X 

4® F(e*>  d5  • V2<x) 

fu£F(^)d4 

L  0 

d 

0 

(4.80) 


Equations  (4.60)  and  (4.61)  give  Vj  and  v2  in  terms  of  u}  and  u2.  Since  the  integrals  in 


equations  (4.79)  and  (4.80)  both  vanish  when  x  =  0,  initially  one  has 


up(0)  =  0 


a 


dx 


x=0 


(4.81) 


Equations  (4.79)  through  (4.81)  represent  the  normal  or  Floquet  form  of  the  first-order 
particular  solution,  its  particular  growth  rate,  and  the  particular  initial  conditions.  These 
can  be  combined  with  the  Floquet  form  of  the  homogeneous  solutions  to  write  a  complete 
integral  for  the  first-order  equation  in  normal  Floquet  form.  At  the  upper  limit  of  the 

integrals  in  equation  (4.79)  when  4  =  x,  the  two  exponential  terms  contained  in  the  two 
Floquet  variables  when  they  are  multiplied  together  become  unity  since  y  and  y2  have  the 

same  magnitude  but  opposite  signs  for  the  case  a  >  1.  If  the  corresponding  <J>.(x),  uQ(x) 

and  F(x,e)  are  bounded  in  [0,  x],  then  as  x  ->««»,  the  first  integral  in  (4.79)  also  approach¬ 
es  °°  since  yl  >  0  and  the  exponential  part  of  v}  dominates  the  solution.  In  the  same 

regard,  the  second  integral  approaches  a  constant  because  y2  <  0  and  the  exponential  part 
of  v2  vanishes.  This  means  that  for  values  of  a  >  1,  as  is  the  case  for  this  analysis,  the 

first  integral  in  equation  (4.79)  is  unstable  while  the  second  integral  is  stable  as  x 

Exactly  how  unstable  the  first  integral  will  become  depends  upon  the  magnitudes  of  die 
parameters  such  as  y}. 
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4.4  First-Order  Equation  Solution 

With  the  homogeneous  solution  and  particular  integral  for  the  first-order  equation 
found  in  both  numerical  and  normal  Floquet  form,  they  can  be  combined  to  form  the  com¬ 
plete  solution.  One  does  this  by  multiplying  the  two  linearly  independent  fundamental 
functions  by  arbitrary  constants  and  adding  these  to  the  particular  solution.  Thus,  the 
complete  solution  to  the  normalized  first-order  equation  in  (4.4)  has  the  form, 
u(x)  =  Auj(x)  +  Bu2(x) 


X 

X 

-Uj(x) 

^F<e,5)d5  +  u2(x) 

«§  F(E-y « 

( 

_ 

(4.82) 


where  A  and  B  are  the  arbitrary  constants  that  must  be  evaluated  from  the  initial  condi¬ 
tions.  It  is  noted  that  although  the  variables  to  the  left  of  the  equality  Uj  and  u2,  they  form 
the  fundamental  solution  set  for  the  first-order  equation,  they  should  not  be  confused  with 
the  notation  in  section  2.4  above  for  the  two  scale  equations  of  the  original  Rayleigh- 
Plesset  equation. 

The  initial  conditions  for  the  first-order  equation  woe  derived  from  the  flaccid 
bubble  analysis  as  in  equation  (2.45)  giving  the  initial  growth  velocity.  From  equations 
(4.5)  and  (4.6)  above  we  have 
Ui(0)  =  0 

,  .  „  -  *W'(1  +  (-  d§)cp=.K  5  C 

It  was  shown  in  die  particular-solution  analysis  that  both  the  particular  solution 
and  its  derivative  are  equal  to  zero  when  x  =  0.  Applying  this  fact  and  equation  (4.83) 
to  equation  (4.82)  and  its  first  derivative,  one  finds  the  equations, 


dui 

diT 


AujtO)  +  Bu2(0)  =  0, 


(4.84) 


AU](0)  +  Bu2(0)  =  C. 


(4.85) 
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From  the  initial  conditions  of  fundamental  set  of  equations  (4.22),  we  And  that  these  two 
equations  yield 


A  =  0, 


B  =  C  . 


(4.86) 


(4.87) 


Substituting  (4.86)  and  (4.87)  into  equation  (4.82),  one  finds  that  the  complete  first-order 
solution  looks  like 


f  n?(£)  f  ui(£) 

u(x)  =  Cu2(x)  -  >.2  Uj(x)  F(e,4)  d^  -  u2(x)  ^^F(e£)dS  ,  (4.88) 


where  C  is  given  by  equation  (4.83).  Differentiating  equation  (4.88),  one  also  finds  that 
the  complete  first-order  growth  rate  is 


i,(x)J^F(e,?)d?-i2(x)  f^|F(e,5)d5  .  (4.89) 


To  express  the  complete  solution  and  its  derivative  in  normal  Floquet  form,  the  u/s 
must  be  expressed  in  terms  of  the  Floquet  variables.  The  substitution  of  equations  (4.77) 
and  (4.78)  into  (4.88)  and  (4.89)  yields 

u(x)  =  C[(X!  -  an) V](x)  +  (X2  -  an)v2(x)] 


X  X 

+  X2a12(X,  -  te)  v,(x)  ^|jF(e,5)d§-v2(x)  ^|jF(e,^)d5  ,  (4.90) 

6  5 


and  its  derivative  is 
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-  C  [(X!  -  a^vrfx)  +  (Xj  -  an)v2(x)] 


+  X2a12(X1-X2)  V!(x)  ^|f(e£)  d£- v2(x)  ^|F(e,4)d4  ,  (4.91) 

L  o  o  J 

where  Vj(x)  =  eYl*  <J>!(x)  and  v2(x)  =  e^x  f2(x)  as  noted  previously  in  equation 

(4.51).  These  are  the  complete  first-order  solution  and  its  derivative  in  terms  of  the  normal 
Floquet  variables,  vj. 

To  study  the  complete  first-order  solution  it  is  best  to  use  the  normal  or  Floquet 
form  which  can  be  continued  from  the  first  period  for  many  periods  in  x.  A  computer 
code  was  developed  which  uses  the  zero-order  solution  and  the  forcing  function  in  terms 
of  the  normalized  bubble  time  x,  along  with  values  of  uv  as  input 

The  first  step  needed  is  the  numerical  calculation  of  the  homogeneous  solutions, 
uj(x)  and  u2(x),  for  one  period  as  0  &  x  £  1.  As  stated  earlier,  this  was  done  by  first 
finding  Zj  and  z2  and  converting  them  to  Uj.  It  was  noted  that  equation  (4.1 1)  becomes 
stiff  as  portions  of  p(x)  become  large  with  increasing  values  of  Uy.  For  cases  of  large  uv, 
standard  numerical  techniques  such  as  a  Runge-Kutta  method  prove  to  be  insufficient  But 
according  to  Gear  [10],  the  fourth-order  explicit  Runge-Kutta  method  produces  convergent 
and  stable  solutions  for  real  values  of  the  quantity  h*p(x)  between  0  and  approximately 
2.7.  Here  h  is  the  step  size  in  x  and  p(x)  comes  from  equation  (4.1 1).  From  equation 


(4.12)  one  sees  that 


so  that 


p(0)  =  2X2  £  p(x)  , 


h*p(x)  £  2hX2  <  2.7  . 


(4.92) 


(4.93) 
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If  200  steps  are  taken  over  the  range  0  £  x  £  1,  equation  (4.93)  yields  conservative  values 
of  X  <  16.4  which  corresponds  to  a  largest  value  of  Uv  «  2.14.  If  100  steps  are  taken  a 
value  of  X  <  1 1.6  is  obtained,  corresponding  to  a  value  of  Uy  *  1.77.  Using  a  variable 

step  size  fourth-order  Runge-Kutta  method,  one  sees  that  the  number  of  steps  taken  over 
the  period  of  x  for  the  entire  range  of  uv  fell  between  100  and  200.  This  then  tells  us  that 
the  numerical  solutions  obtained  from  this  method  are  stable  for  values  of  uv  at  least  up  to 
1.77. 

By  looking  at  the  curves  of  p(x)  in  figure  4.1  for  higher  values  of  uv,  for  a  number 
of  steps  between  100  and  200,  we  see  that  the  value  of  p(x)  only  violates  the  inequality  in 
equation  (4.93)  over  very  narrow  strips  near  the  values  of  x  =  0  and  x  =  1.  These  regions 
of  stiffness  produce  unstable  solutions  over  very  small  intervals  compared  to  the  rest  of  the 
range  of  x  where  stable  solutions  are  found.  It  was  felt  therefore,  that  the  errors  accumu¬ 
lated  using  this  variable  step  size  Runge-Kutta  method  might  be  acceptable.  This  is  indeed 
the  case,  which  will  be  verified  later  when  the  Wronskian  of  the  system  is  calculated  for  a 
range  of  uv  values.  Thus,  we  are  quite  fortunate  in  the  fact  that  even  hough  there  exist 
regions  of  numerical  instability  over  the  range  of  x  for  the  solutions  to  equation  (4. 11), 
these  regions  are  small  and  a  reliable  proven  numerical  technique  can  be  used  to  find  the 
solutions  with  acceptable  error. 

The  use  of  the  WKB  method  was  explored  for  this  problem  because  the  value  of 
the  parameter  A.  becomes  large  as  u*  increases.  It  was  found  that  although  the  equation  is 

certainly  stiff  in  such  cases,  many  terms  beyond  the  order  of  1  A2  must  be  obtained  if  one 
is  to  get  a  sufficiently  accurate  asymptotic  expansion. 

The  variable  step  size  fourth-order  Runge-Kutta  method  was  applied  once  for  each 
solution,  z1  and  Z2,  and  these  were  then  converted  back  to  u/s  using  equation  (4.10). 
These  solutions  and  their  derivatives,  for  one  value  of  Uy,  are  shown  in  figures  4.3. 

Later  on  in  the  analysis  the  values  of  such  solutions,  along  with  uq(x),  will  need  to  be 
given  in  even  increments  of  x  for  numerical  integration  purposes.  To  obtain  evenly 
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spaced  values,  a  cubic  interpolation  scheme  was  used,  which  picked  a  target  value  of  x 
and  used  the  surrounding  four  known  values  of  the  solution  to  interpolate  the  corre¬ 
sponding  value.  The  use  of  an  interpolation  scheme  is  also  apt  to  add  error  to  the 
solution.  But  since  the  solutions  are  well  behaved  and  a  third-order  method  [15]  was 
used,  these  errors  were  felt  to  be  insignificant.  This  same  method  was  applied  to  get 
evenly  spaced  values  of  Uq  as  input  to  the  calculation  of  the  first  order  solution. 

The  values  of  the  fundamental  solutions  and  their  derivatives  at  x  =  1  are  next  used 
to  determine  the  constants  in  the  matrix  [A]  from  equation  (4.26).  From  these  values  the 
Wronskian  of  the  system  can  be  determined  using  equation  (4.35).  This  Wronskian  is 
used  to  assess  the  accuracy  of  the  numerical  calculations  of  the  fundamental  solutions.  As 
was  shown  for  Hill's  equation,  the  Wronskian  should  have  a  constant  value  of  unity.  By 
monitoring  the  value  of  the  Wronskian  for  given  values  of  Uy,  the  confidence  that  the 

solution  is  sufficiently  accurate  is  measured.  Figure  4.4  shows  the  scatter  in  values  of  the 
error  (A  - 1)  over  the  range  of  Uy.  It  is  seen  that  although  the  variations  from  A  »  1 
increase  with  Uy,  the  difference  is  less  than  10*5,  which  is  acceptable  for  this  analysis. 
Closer  inspection  of  this  figure  shows  that  great  accuracy  is  obtained  for  values  of  Uy  up  to 
around  2,  which  verifies  the  numerical  stability  of  the  u,  solutions  talked  about  earlier 
using  the  work  of  Gear  [10].  As  the  value  of  Uy  gets  larger,  the  stiff  regions  of  the 
equation  produce  some  error  in  the  solutions,  but  not  enough  to  discredit  the  present  use  of 
the  variable  step  size  Runge-Kutta  method. 

The  stability  parameter,  a,  can  also  be  determined  from  the  values  of  matrix  [A] 
using  equation  (4.34).  The  value  of  this  parameter  leads  to  dynamic  stability  assesment  of 
the  fundamental  set  of  solutions  used  in  the  complete  first-order  solution.  The  fundamental 
set  of  solutions  begins  at  a  neutral  stability  point  when  Uy  =  1,  corresponding  to  a  value  of 

a  =  1,  as  shown  in  figure  4.5.  This  figure  also  shows  that,  exept  for  the  special  value  of 
uv=  l.ais  greater  than  unity  over  the  whole  range  of  Uy  values.  Thus,  for  1  >  Uy  >  1, 

a  >  1  and  the  fundamental  set  of  solutions  has  one  bounded  and  one  unbounded  solution 
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Stability  parameter,  a. 
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Figure  4.5  Stability  of  the  First-Order  Fundamental  Set  of  Solutions  for  a  Range  of  uv 
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as  rime  increases,  and  excupt  for,  uv  =  1,  no  values  of  Uy  do  both  of  the  the  solutions 
become  stable.  For  the  case  of  uv  =  1,  a  =  1  and  ^  =  ^2  =  1,  the  fundamental  set  reduces 
to  neutrally  stable  trigonometric  solutions  ofiinit  period  and  becomes  unstable  again  with 
further  decreases  in  uv. 

The  next  step  is  to  determine  the  Floquet  variables,  v;(x)  and  Vj(x).  For  computa¬ 
tional  puiposes,  it  is  desired  have  them  in  the  form  of  equation  (4.51)  which  involves  a 
well  defined  exponential  function  multiplied  by  a  periodic  function.  From  the  a  value,  X] 
and  X2  can  be  determined  from  equation  (4.41).  From  these,  Yi  and  Y2  can  be  determined 
from  equation  (4.41)  with  a  unit  period  of  T  =  1.  The  known  characteristic  exponents  take 
care  of  the  exponential  part  of  the  Floquet  variables. 

Values  of  v;  for  one  period  in  x  can  be  found  using  the  fundamental  set  of 
solutions  and  equations  (4.60)  and  (4.61).  In  the  same  fashion,  v;  can  be  found.  Then 

using  equation  (4.62)  and  its  derivative  with  respect  to  x,  the  periodic  <J>j  and  <j>;  can  be 
determined.  Figure  4.6  shows  a  plot  of  and  for  several  values  of  Uy.  This  then 
takes  care  of  the  periodic  part  of  the  Floquet  variables.  There  now  exists  a  form  of  v;(x) 
and  Vj(x)  that  can  be  continued  over  many  periods  in  x. 

It  remains  now  to  determine  the  integrals  involving  the  forcing  function  from  the 
particular  solution.  This  is  simply  done  using  Simpson's  rule  for  evenly  spaced  values  of 
u0(x),  V;(x),  and  F(ept).  The  same  process  is  applied  to  the  integrals  in  the  first  derivative 

equation.  All  these  components  are  then  combined  as  in  equations  (4.90)  and  (4.91).  The 
constant  C  can  be  determined  using  the  forcing  function  and  equation  (4.83). 

The  above  steps  are  combined  into  a  computer  program  contained  in  appendix  A. 
Some  of  the  input  needed  for  the  code  is  the  pressure  coefficient  data  for  a  circular  cylinder 
and  a  cavitation  number.  For  this  example,  a  value  of  K  =  2.3  was  used  along  with  a 
typical  value  of  the  free  stream  velocity  of  40  fps.  The  initial  nucleus  radius  size  was  taken 
to  be  7  microns  and  a  characteristic  length  for  the  circular  cylinder  was  D  =  2  in.  Values  of 
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the  surface  tension  and  density  of  water  at  70  degrees  were  used  and  a  dissolved  air 
saturation  pressure  of  10  psia  was  also  used.  The  computed  value  of  uv  for  these  input 
parameters  is  Uy  =  0.5994  and  e  is  l.tixlO-6.  But  figures  4.7  and  4.8  show  the  complete 
first-order  solution  and  its  derivative  over  many  periods  in  x  for  a  value  of  uv  =  1.5994 . 

The  reason  for  this  change  is  that  the  calculated  value  of  Uy  from  these  input  data 
fell  below  the  critical  value  of  unity  for  this  investigation.  The  calculated  value  of  uv  corre¬ 
sponds  to  small  scale  compressive  air  bubble  oscillations  which  characterize  the  zero-order 
phase  plane  trajectories  to  the  left  of  the  Uq(x)  =  1  point  shown  in  figure  3.1.  As  was 

stated  in  chapter  3,  the  present  numerical  analysis  deals  only  with  vaporous  growth 
trajectories  to  the  right  of  this  point.  The  extension  of  the  code  to  smaller  values  of  uv 
could  have  been  made,  but  it  was  not  carried  out;  because  the  dynamic  stability  analyses 
for  values  of  uv  to  the  left  as  well  as  to  the  right  of  the  critical  point  show  that,  aside  from 
the  case  of  neutral  stability  at  Uy  =  1,  all  first  order  solutions  are  unstable  as  shown  in 
figures  4.2  and  4.5.  We  concluded  therefore  that  calculations  carried  out  for  uv  =  1.5994 
can  give  a  useful  illustration  of  the  consequences  of  the  first-order  instability. 

It  is  seen  from  figure  4.7  that  the  total  ut  solution  grows  exponentially  with  x.  For 
this  case  the  value  of  a  is  greater  than  1,  thus  producing  exponential  terms  with  positive 
exponents  that  will  become  unbounded  as  x  increases.  Figure  4.8  shows  the  derivative  of 
this  complete  integral.  Figures  4.9  and  4.10  show  the  homogeneous  and  particular  parts  of 
the  first-order  solution.  From  them  we  see  that  the  large  magnitude  of  the  total  Uj  solution 

comes  from  the  homogeneous  part.  This  is  due  the  constant  C  multiplying  the  comple¬ 
mentary  function.  C  is  of  the  order  103  for  the  given  input  values.  Because  this  constant 
term  is  relatively  large,  it  dominates  the  smaller  oscillations  from  the  particular  part  of  the 
solution.  The  particular  part  of  the  solution  grows  at  a  slower  rate  than  the  homogeneous 
part  due  to  the  fact  that  die  amplitude  of  the  forcing  function  contained  in  the  particular 
integrals  grows  rather  slowly  in  terms  of  bubble  time  so  that  the  modulated  first  order 


iXi 


Figure  4. 10  Particular  Part  of  the  First-Order  Solution  as  a  Function  of  Normalized 
Bubble  Time 
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oscillations  show  a  less  rapid  growth  than  does  the  complementary  solution.  But  because 
ofthe  presence  of  the  complementary  function  in  the  particular  solution,  it  is  also  basically 
unstable.  It  must  be  remembered  also  that  since  e  is  of  the  order  of  10*6 ,  the  entire  first 
order  solution  will  be  of  the  order  of  10"3  compared  to  the  zero-order  solution.  But  it 
appears  from  their  forms  that  the  present  zero  and  first  order  solutions  can  not  lead  to  any 
useful  inception  criteria.  It  is  necessary,  therefore,  that  the  suggestion  developed  at  the  end 
of  section  3.53  above  be  seriously  considered  in  a  separate  investigation.  This  is  particu¬ 
larly  important  because  it  appears  that  the  zero-order  solution  will  frequently  dominate  the 
first  order  solution,  at  least  in  the  initial  phases  of  the  forced  bubble  motions  as  examined 
here.  Therefore,  a  more  general  reexamination  of  the  zero-order  solution  centered  around 
the  integrability  conditions  is  warranted. 

The  preceding  results  do  not  restrict  the  value  of  Uy  and  for  that  reason  it  has  not 
been  possible  to  give  a  closed  form  perturbation  solution.  On  the  other  hand  when  uv  =  1 
the  solution  is  neutrally  stable  and  there  is  a  possibility  that  the  farcing  function  can  have 
more  influence  cm  the  bubble  growth  than  was  found  in  die  proceeding  example.  Therefore 
it  is  of  interest  to  obtain  a  solution  for  this  special  case  even  though  the  condition  uv  =  1  is 
too  severe  a  restriction  to  make  the  result  of  practical  value. 

4.4.1  The  Solution  when  u„  =  1 

In  Section  3.4.1  we  found  in  the  limit  as  uv  -» 1,  that  die  zero-order  oscillation 
parameter  is  X  =  JtV2  and  according  to  Equation  (2.20)  the  period  is  T  =  (l+q)3/4A. . 
Therefore  when  uv  «  1, 1 + q  =  y  and  T  =  7ty3/4V2.  Moreover,  the  zero-order  au¬ 
tonomous  differential  equation  (3.1)  and  its  initial  conditions,  (3.2)  and  (3.3),  insure  that 
at  uv  *  1  the  zero-order  solution  is  simply  uq(t)  =  1  for  all  values  of  the  bubble  time  x,  as 
discussed  in  connection  with  Equation  (3.14).  Therefore  we  are  free  to  consider  the  first 
order  differential  equation  (4.4),  for  Uj(x),  and  its  initial  conditions,  equations  (4.5)  and 
(4.6).  The  fact  that  Uq  =  constant  in  this  case  causes  die  derivatives  of  uo  to  vanish  in 
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equation  (4.4)  and  the  limiting  value  of  X  together  with  the  unit  value  of  Uq  permits  one  to 
simplify  equation  (4.4)  to  read 

^  +  4tc2Ui  =  2rc2F(ts).  (4.92) 

In  equation  (4.92)  the  forcing  function  argument  has  been  transformed  from  real  time,  t,  to 
dimensionless  laboratory  time  tg  =  tD/V0 ,  so  that  now 

F0s)=  4  Wr  [-  Cp(tj)  -  K]  =iwrF,(y 

and  in  accordance  with  the  discussion  of  Appendix  B,  the  term  Fjftg)  represents  a  closed 
form  approximation  to  the  numerical  data  used  above.  In  shifted  laboratory  time  Fj  is  given 
by  equation  (B.3).  Moreover,  the  normalized  bubble  time,  x,  is  now  a  reduced  form  of 
equation  (3.17)  which  becomes 

X=7cy3/4V2  *  <4-93) 

Since  tg  =  er,  it  follows  that  the  laboratory  time  and  the  normalized  bubble  time  are  related 
by 

tj  =  e  x  n  y3/4V2  .  (4.94) 

Consequently  when  x  is  the  independent  variable  the  initial  condition  of  equation  (4.83) 
above  becomes 

(495) 

while  we  have  the  condition  Uj(0)  =  0  as  before. 

The  solution  of  Equation  (4.92),  subject  to  the  prescribed  initial  conditions  is 
well  known: 

Uj  =  ^  sin  2joc  +  j  WrJ^(e;  §)  sin  2x(x-^)  d£ ,  for  x  £  xs,  (4.96a) 


and  if  x  >  x,, 
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a8 

U!  sin2jtx  +  jWr‘  J^(e;£)sin2jc(x-5)dl; 

. 

+  (2.225  -  K)fsin27t(x-5)dS  >.  (4.966) 

xs 

If  we  let  b  =  7ty3/4V2  we  can  write  the  approximate  formula  for  Fj(e;  x)  for 
tj  ^  ts  or  for  x  £  xs,  as 

Fj (e;  x)  =  aj-K  +  a2to(K)  +  a2ebx + a3sin[a4ebx  +  a4to(K)+  a5] ;  (4.97) 

and  in  the  laminar  separation  zone  when  x  >  xs  =  (0.365  -  to)/(eb), 

Fjfc;  x)  =  2.225  -  K,  (4.98) 

as  already  indicated  in  equation  (4.966 ).  The  quantity  to(K)  is  the  origin  of  the  forcing 
function  and  is  given  analytically  in  figure  B4.  The  quantities  ?s  and  xs  denote  the 
separation  point  in  shifted  coordinates,  t  =  ts-  Xq. 

Because  of  the  form  of  equation  (4.97),  the  integral  in  equation  (4.96a)  can  be 
expressed  as  the  sum,  Ij  + 12  + 13,  for  'i  £  or  x  <.  xs.  If  x  is  greater  than  xs  there 
also  will  be  a  fourth  integral,  I4,  given  by  the  second  integral  of  (4.966).  Considering 
these  integrals  in  order,  one  finds  that 

Ij  =  — K  2tc2^K~ °°s 27CX)’  x  ^  xs  or 

Ij  =  — K  2/^-  [cos2tc(x-xs)  -  cos2tcxs],  if  x  >  Xs. 

a2cb  ri  1 

h=  2ir  [x’2Ksin2nx  J.  tfxSxs  or 

I2  =  ^xscos2jt(x-xs)  +  2^(sin27t(x-xs)-sin2jtx]  J,  if  x  >  xs. 

Using  the  trigonometrical  addition  formula;  and  integrating,  one  can  write  I3  as 
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i  -ii{  — L 
3  2  [2jc+a4eb 


[  sin(a4ebx  +  <j>)  +  sin(2jtx  -<(>)] 


1 


[sin(a4£bx+<|>)  -  sin(2jtx  +  <)>)]  l , 


2 n - a4eb 

where  <)>  =  a4to(K)  +  a5  and  0  <  x  <  xs.  For  the  case  in  which  x  £  xs  we  have, 


1 


2  [2jt  +  a4eb 


[sin{a4£bx  +  <|>  -27t(x-Xs)}  +  sin (2tcx -<)))] 


2?c 


[sin{a4£bx  +  <{>  +  27t(x-Xs)}  -  sin(2rcx  +  <())]| . 


Continuing  with  the  case  x  >  xs,  we  evaluate  the  last  integral,  I4  in  equation 
(4.96b ),  in  order  to  account  for  the  constant  pressure  in  the  separation  bubble.  This  last 
result  is 


T  2.225 -K  r,  „  , 

I4  =  — 2tc — [l-cos2*(x-Xg)]. 

Then  the  complete  closed-form  solution  to  order  £  can  be  written  as 

u(x)  =  Uq  +  £U!  =  1  +  £  |  ^  sin2nx  +  J  Wr [  Ij  + 12  + 13+ 14]  J,  (4.99) 

where  I4  =  0  if  x  £  xs .  Because  of  the  approximations  representing  Fj(£,x)  from 

Appendix  B  and  the  stability  of  the  solution  in  this  case,  one  can  examine  the  influence 

of  the  forcing  function  on  the  solution  given  in  equation  (4.99).  As  already  found  in  the 

general  case,  the  term  from  the  initial  condition  simply  provides  a  steady  oscillation  of 
eC 

amplitude—  at  the  natural  frequency  of  the  bubble. 

The  interesting  part  of  the  solution  is  found  in  the  sum  of  the  integrals.  We  shall 
use  the  artificial  example  of  e  =  0.001,  K  =  2.3,  y  =  12.9  and  b  =  30.3.  From  these  data 
and  with  the  help  of  equation  (4.94)  and  figure  B5  one  sees  that  xs  =  10.  We  also  find  that 
to  =  0.06  from  figure  B4.  Then  since  we  have  y,  K  and  uv,  the  Weber  number  can  be 

found  from  figure  3.2,  and  it  is  Wr  =  24.  These  parameters  and  the  coefficients,  al, . . .,  a5, 
from  figure  B3  are  all  that  one  needs  in  order  to  evaluate  the  term,  ^Wr[  Ij+ 12+ 13+  I4],  in 

equation  (4.99).  The  result  of  this  calculation  is  shown  in  figure  4.1 1,  below.  The  figure 
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shows  that  the  fundamental  vibration  of  the  bubble  is  modulated  by  the  forcing  function  as 
one  expects.  The  bubble  grows  and  starts  to  collapse  but  when  the  bubble  reaches  the 
laminar  separation  zone  its  collapse  is  anrested  and  the  bubble  simply  vibrates  in  place.  The 
growth  shown  by  this  example  is  modest,  compared  to  the  initial  condition  which  con¬ 
tributes  the  term,  43.35  sin2jcx.  Consequently,  even  though  this  initial  term  provides  no 
growth  its  effect  will  obscure  the  forcing  function  response.  It  also  evident  that  the 
condition  uv  =  1  is  too  restrictive.  No  inception  calculations  can  be  earned  out  if  uv  is  so 
constrained  even  if  the  growth  were  greater.  As  this  example  shows,  however,  the  first 
order  solution  does  respond  to  the  forcing  function  as  expected. 


Figure  4.1 1  Particular  First-Order  Solution  for  uv  =  1  as  Calculated  from 
the  Convolution  Integrals  of  Equations  (4.96a)  and  (4.96 b). 
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CHAPTER  5 

SUMMARY  AND  CONCLUSIONS 


This  investigation  has  sought  to  solve  the  forced  Rayleigh-Plesset  equation  for  the 
dynamics  of  forced  vaporous  growth  of  a  cavitation  bubble  in  an  inviscid  incompressible 
flow  about  submerged  bodies  using  a  two-scale  perturbation  expansion.  Of  course,  the 
introduction  of  a  viscous  term  into  the  equation  could  well  introduce  a  third  time  scale. 
Experience  with  straight  forward  numerical  integrations  of  the  forced  Rayleigh-Plesset 
equation  with  and  without  viscosity  has  shown  the  viscous  term  to  have  negligible  effect 
for  relatively  long  pulse  forcing  functions  of  interest  here.  Other  simplifying  assumptions 
were  that  the  isothermal  bubbles  are  spherical  and  that  they  do  not  interact  with  each  other 
or  with  the  surface  of  the  submerged  body.  For  further  simplification  we  suppose  that  the 
pressure  distribution  on  the  body  supplies  the  excitation  for  the  bubble  dynamical  equation 
and  that  the  boundary  layer  slows  the  progress  of  the  bubble  along  the  surface  in  a  very 
rudimentary  way.  The  bubble  growth  is  not  coupled  to  the  shear  layer.  The  two  time  scales 
arise  firstly,  from  the  natural  period  of  a  microbubble  itself  and  secondly,  from  the  time 
required  for  the  growing  vapor  bubble  to  move  through  the  low  pressure  region  on  the 
body.  The  first  of  these  dimensionless  times  is  the  fast  "bubble  time"  and  the  second  is  the 
slow  "laboratory  time".  The  ratio  of  slow  to  fast  times  is  the  perturbation  parameter,  e. 

The  first  part  of  this  investigation  dealt  with  the  forcing  function  that  drives  the 
growth  of  the  bubble.  Baker  [2]  used  pressure  coefficient  data  from  a  hemispherical 
hcadform  having  a  short  laminar  bubble  followed  by  turbulent  reattachment.  The  present 
analysis  considers  a  new  possibility  for  a  forcing  function,  namely  the  supercritical  flow 
about  a  circular  cylinder,  because  this  flow  exhibits  a  laminar  separation  bubble  just 
upstream  of  the  turbulent  separation  point  It  presents  an  instance  of  cavitation  inception 
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not  previously  studied  in  this  context,  and  may  help  to  determine  if  this  type  of  analysis  can 
be  applied  to  any  body  which  has  a  boundary  layer  with  a  short  laminar  separation  bubble. 
Using  pressure  coefficient  data  for  a  circular  cylinder  and  employing  a  very  simple  bubble 
convection  model ,  we  converted  these  data  into  a  forcing  function  dependent  on  the  slow 
laboratory  time. 

Next,  the  kinematics  of  a  bubble  as  it  travels  through  the  high  pressure  region  on 
the  body  were  considered  because  this  region  lies  directly  upstream  of  the  vaporous  growth 
region.  These  results  were  used  to  derive  initial  conditions  at  the  start  of  vaporous  growth. 
Using  the  two  time  scales  of  the  present  problem,  the  dynamical  equations  and  their  initial 
conditions  were  derived  to  the  third  order  of  the  small  parameter,  e.  A  few  errors  were 
found  and  corrected  in  similar  equations  derived  by  Baker. 

The  relationships  between  the  measurable  flow,  body  and  bubble  parameters  were 
then  considered  in  order  to  indicate  their  effect  on  the  range  £.  It  was  found  useful  to 
express  these  several  parameters  as  dimensionless  quantities  and  to  relate  them  to  practical 
experimental  ranges  of  the  variables  to  cover  likely  applications  of  this  analysis.  As  Baker 
showed  from  his  zero-order  phase  plane  study,  there  is  a  single  parameter,  i^,  the  affix  of 

the  vortex  point  in  the  zero  order  phase  plane,  that  relates  the  bubble  dynamics  to  key  flow 
parameters  such  as  the  cavitation  number,  Weber  number  and  dissolved  air  content  This 
then  is  the  single  parameter  which  characterizes  each  member  of  the  family  of  solutions. 

A  solution  of  the  zero-order  equation  was  next  obtained.  For  small  values  of  uv  in 

the  neighborhood  of  unity,  the  zero-order  solution  is  nearly  sinusoidal.  Thus  a  simple 
closed  form  approximate  solution  for  small  oscillations  in  Uy  was  found  using 
trigonometric  functions  which  correctly  displayed  the  periods  and  amplitudes  of  the  u0 
solution.  Although  Baker’s  zero-order  solution  in  the  form,  x  =  f(uo),  showed  good 
results,  it  can  not  be  inverted,  The  solution  of  the  first-order  equation  requires  the  zero 
order  solution  to  have  the  form  uo(x).  Therefore,  a  numerical  uo  solutions  were  found  for 
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one  period  in  x  that  arc  accurate  over  the  range  of  uv  values.  The  solution  for  one  period 
can  be  reused  over  large  ranges  in  x  due  to  its  periodicity. 

Considerable  study  centered  upon  the  first-order  integrability  condition  and  it 
suggests  that  the  zero-order  solution  should  be  independent  of  the  slow  laboratory  time. 
The  integrability  condition  is  a  second  order  partial  differential  equation  for  u0  consisting  of 
those  terms  having  mixed  partial  derivatives  in  the  first  order  equation  that  must  be  set 
equal  to  zero.  The  above  independence  will  be  true  only  if  the  arbitrary  functions  of  the 
slow  time,  which  result  from  partial  integrations  with  respect  to  the  fast  time,  are  in  fact 
completely  independent  functions.  This  is  the  situation  ordinarily  encountered  and  it  is  a 
fundamental  precondition  in  thpresent  investigation.  That  the  zero  order  solution  is  then 
independent  of  ^  was  shown  explicitly  for  the  small  oscillation  approximation  to  the  zero 
order  solution.  Next  the  integrability  condition  itself  was  examined  from  the  viewpoint  of 
the  theory  of  charactersitics.  In  addition  the  most  general  functional  form  of  u0  was  found 
in  terms  of  the  two  time  scales  by  direct  integration.  These  considerations  verified  that  the 
zero  order  solution  should  not  depend  on  independent  arbitrary  functions  of  the  slow 
time.  The  importance  of  this  comes  in  the  fact  that  all  terms  involving  derivative  of  uq  with 
respect  to  t,.  in  the  higher  order  equations  vanish,  presumably  suppressing  secular  terms  in 
the  solutions.  It  also  allows  both  the  zero  and  higher  order  order  equations  to  be  written  as 
ordinary  differential  equations  instead  of  partial  differential  equations. 

Finally,  the  first-order  equation  was  solved.  This  turned  out  to  be  a  nonhomo- 
geneous  linear  second-order  differential  equation  with  periodic  coefficients  which  are 
functions  Uq.  It  was  found  that  the  homogeneous  form  of  this  equation  could  be 
transformed  into  a  special  form  of  Hill's  equation.  For  small  oscillations  in  Uy,  the  first 
order  solution  can  be  approximated  in  terms  of  Mathieu  functions.  This  fact  was  used  to 
investigate  the  departure  of  the  solution  from  a  neutrally  stable  state,  involving  only  simple 
harmonic  oscillations,  to  unstable  states  as  the  oscillation  amplitude  increases  with 
increasing  Uy. 
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The  linearity  of  the  first-order  equation  enabled  the  superposition  of  the  compli¬ 
mentary  solution  and  the  particular  integral  in  order  to  arrive  at  the  complete  integral  for  the 
first-order  problem.  The  solution  to  the  homogeneous  form  of  the  first-order  equation  can 
be  found  using  Floquet  theory.  This  solution  in  normal  Floquet  form  was  found 
numerically  for  one  period  in  x,  then  extended  over  more  periods  using  the  periodicity  of 
the  Floquet  variables.  Variation  of  parameters  was  then  used  to  derive  the  particular 
integral  and  the  initial  conditions  were  used  in  conjunction  with  the  complete  solution  to 
evaluate  the  constants  of  integration.  Again  it  was  found  that  these  constants  are 
independent  of  laboratory  time,  and  that  only  the  bubble  time  appears  explicitly  in  the 
solution. 

Normally  the  Uj  solution  would  be  multiplied  by  the  first-order  of  e  and  be  added  to 
the  uo  solution  to  obtain  an  approximate  result  to  the  governing  Rayleigh-Plesset  equation. 

But  in  this  investigation,  only  the  nature  of  the  complete  first-order  solution  was  obtained. 
This  is  because  the  actual  value  of  Uv  which  is  consistent  with  the  other  flow  parameters  for 
the  test  case  is  =  0.5994.  This  value  defines  a  zero  order  phase  plane  trajectory  to  the 
left  of  Uv  =  1.0,  denoting  compressive  oscillations  of  Uq.  The  present  study  has  been 
limited  to  cases  for  which  the  zero-order  trajectories  lie  to  the  right  of  uv  =  1.0, 
corresponding  to  oscillations  of  expansion  which  would  naturally  seem  to  lead  to  further 
vaporous  growth  under  the  influence  of  a  forcing  function.  The  present  example  shows 
that  this  may  not  always  be  true.  Nevertheless  a  specific  example  for  uv  =  1.5994  has  been 
worked  out  in  order  to  indicate  important  features  of  the  present  solution. 

New  things  that  have  emerged  from  this  work  include  the  small  oscillation  theory 
which  showed  that  trigonometric  functions  can  be  used  to  derive  approximate  results  to 
both  the  zero  and  First-order  solutions  for  values  of  U*  close  to  unity.  To  the  authors' 

knowledge,  this  is  the  first  application  of  the  Floquet  theory  to  bubble  dynamics.  This 
approach  makes  it  possible  to  investigate  the  dynamic  stability  of  the  solutions  and  we 
found  that  aside  from  die  small  oscillation  case  at  Uy  =  1,  all  other  solutions  are  unstable.  It 
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was  also  useful  because  numerical  work  was  only  needed  over  the  first  period  of  the  first 
order  complementary  complementary  function,  Uj,  thereby  minimizing  error  accumulation 

since  the  complete  solution  is  extended  over  many  periods.  The  forcing  function  for  the 
circular  cylinder  was  used  to  show  one  example  of  the  solution.  It  is  presented  to  demon¬ 
strate  that  a  numerical  solution  can  be  found  using  different  but  appropriate  bodies.  The 
chief  practical  physical  interpretation  to  be  obtained  from  these  calculations  is  that  the  lack 
of  stable  solutions  prevents  their  use  for  the  inception  problem. 

Finally  we  note  that  our  our  decision  to  satisfy  the  first  order  integrability  condition 
by  taking  the  zero  order  solution  to  be  independent  of  the  slow  laboratory  time,  ts,  has  had 
a  profound  effect  on  the  outcome  of  this  research.  This  choice  was  based  upon  the  fact  that 
in  the  absence  of  additional  conditions  on  the  problem,  constants  or  functions  of  integration 
are  arbitrary  independent  functions  of  the  slow  time.  Perhaps  the  most  important  effect  of 
this  requirement  may  be  the  dynamic  instability  found  in  the  first  order  solution.  Conse¬ 
quently  it  appears  that  if  one  abandons  his  insistence  that  that  the  functions  of  integration  in 
the  zero  order  solution  be  independent  arbitrary  functions  of  the  slow  time,  then  further 
progress  might  be  possible  and  a  physically  acceptable  solution  might  be  the  outcome. 
Whether  or  not  a  generalized  approach  to  the  problem  as  suggested  in  section  3.5.3  above 
can  succeed  in  changing  the  present  unsatisfactory  state  of  affairs  remains  to  be 
investigated. 
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APPENDIX  A 
FORTRAN  CODE  LISTING 


This  is  a  listing  of  the  computer  code  used  to  do  all  the  calculations.  This 
includes  calculating  the  forcing  function,  zero-order  solution  and  first-order  solution. 
Since  phases  of  the  code  were  developed  at  different  times  throughout  the  analysis  and 
pieced  together  into  this  final  code,  its  working  order  is  by  no  means  the  most  efficient 
The  important  result  is  that  it  works. 

The  code  is  written  in  standard  FORTRAN  77,  and  was  run  on  both  a 
MACINTOSH  PC  and  a  VAX  1 1-780  system.  Many  of  the  numerical  subroutines  have 
been  adapted  from  Numerical  Recipes  by  Press  et  al.  [15]  which  made  much  of  the 
programming  easier  and  more  efficient 

The  inputs  needed  for  the  code  are  provided  from  three  different  files.  The  first 
file  corresponding  to  unit  10  contains  initial  values  to  perform  the  first-order  integration. 
These  include  the  initial  and  final  values  of  x  for  (me  period,  XI  and  X2,  an  estimated 
step  size  to  be  used  in  the  Runge-Kutta  routine,  H,  the  required  accuracy  used  in  the 
integration,  EPS,  the  minimum  step  size  allowed,  HMIN,  and  the  number  of  first-order 
equations  to  be  solved,  N,  which  is  two  for  this  problem.  The  code  then  reads  the 
initial  values  of  the  fundamental  set  of  solutions  into  the  2  dimensional  array  Z0. 

The  next  input  file  corresponding  to  unit  1 1  contains  the  basic  parameters  of  the 
flow.  These  include  the  cavitation  number,  XKVAL,  the  free  stream  velocity  in  ft/s, 
VNOT,  the  initial  free  stream  nucleus  radius  in  itm,  RNOT_MIC,  the  characteristic 
body  length  in  in.,  DIAM_INCH,  die  density  of  the  liquid  in  slugs/ft3,  RHO,  the 
surface  tension  of  the  liquid  in  lb/ft,  SIGMA,  and  the  dissolved  air  pressure  in  psia, 
PAJNCH.  Unit  12  corresponds  to  the  input  file  containing  die  pressure  coefficient 
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versus  dimensionless  arc  length  data.  It  is  read  by  first  inputing  the  total  number  of  data 
points,  NUM_PTS,  and  then  reading  the  table  for  S  and  CP. 

These  parameters  are  input  with  their  given  physical  dimensional  units.  The 
code  automatically  converts  each  parameter  into  consistent  units  and  eventually  into  the 
proper  nondimensional  form.  It  then  converts  the  input  into  a  corresponding  value  of 
Uy  and  forcing  function,  then  it  solves  the  zero-order  and  first-order  equations.  The 

output  of  the  program  includes  the  total  first-order  solution  into  unit  40,  the  first-order 
stability  parameters  into  units  15  and  16,  and  a  summary  of  the  inputed  and  calculated 
flow  parameters  into  unit  45.  Simple  modifications  can  be  made  to  the  code  in  order  to 
output  any  of  the  intermediate  variables. 
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PRESSURE  DATA  TABULATION  AND  FORCING  FUNCTION  CALCULATION 


B  1.  Data  Tabulation 

This  appendix  tabulates  two  sets  of  pressure  coefficient,  Cp,  versus  dimensionless  arc 
length,  s,  for  a  circular  cylinder  at  supercritical  Reynolds  numbers.  Table  B.l  contains  the 
original  data  as  read  from  a  graph  given  by  Gowen  and  Perkins  [4].  Since  these  data,  as 
read,  did  not  produce  a  smooth  curve,  Lagrangian  interpolation  was  used  to  produce  evenly 
spaced  points  over  the  interesting  range  of  s.  These  points  were  then  smoothed.  They  are 
given  in  in  Table  B.2  and  they  are  the  input  Cp(s)  data  used  in  the  Fortran  code. 

Table  B.l  Original  Data  Table  B.2  Smoothed  Data 


No. 

S 

Cp 

No. 

S 

Cp 

1 

0.0000 

1.0000 

1 

0.000 

1.00000000 

2 

0.0403 

0.9600 

2 

0.064 

0.93858008 

3 

0.0873 

0.9200 

3 

0.128 

0.87090363 

4 

0.3491 

0.5000 

4 

0.192 

0.78682947 

5 

0.6981 

-0.5750 

5 

0.256 

0.68326172 

6 

1-0472 

-1.8313 

6 

0.320 

0.55809864 

7 

1.2217 

-2.2531 

7 

0.384 

0.41404022 

e 

1.2654 

-2.3063 

8 

0.448 

0.25C18646 

9 

1.3090 

-2.3813 

9 

0.512 

0.06653735 

10 

1.3526 

-2.4563 

10 

0.576 

-0.13690708 

11 

1.3963 

-2.5000 

li 

0.640 

-0.36078890 

12 

1.4399 

-2.5531 

12 

0.704 

-0.60684652 

13 

1.4835 

-2.5688 

13 

0.768 

-0.86596418 

14 

1.5272 

-2.5667 

14 

0.832 

-1.11726326 

15 

1.5708 

-2.5583 

15 

0.896 

-1.35007730 

16 

1.6144 

-2.5000 

16 

0.960 

-1.56440630 

17 

1.6581 

-2.4406 

17 

1.024 

-1.76025025 

18 

1.7017 

-2.3500 

18 

1.088 

-1.93760915 

19 

1.7453 

-2.2750 

19 

1.152 

-2.09648302 

20 

1.8326 

-2.2438 

20 

1.216 

-2.22826666 

21 

1.9200 

-2.1625 

21 

1.280 

-2.33585874 

22 

1.344 

-2.43371901 

23 

1.408 

-2.51203038 

24 

1.472 

-2.55500333 

25 

1.536 

-2.55466889 

26 

1.600 

-2.50791555 

27 

1.664 

-2.41927707 

28 

1.728 

-2.32443463 

29 

1.792 

-2.26472235 

30 

1.856 

-2.22700000 

31 

1.920 

-2.22500000 

32 

1.984 

-2.22500000 

33 

2.048 

-2.22500000 

34 

2.112 

-2.22500000 

35 

2.176 

-2.22500000 

36 

2.240 

-2.22500000 

37 

2.304 

-2.22500000 
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B  2.  Pressure  Distribution  and  Forcing  Function  on  the  Cylinder 

As  noted  in  Chapter  2  we  have  used  the  data  as  read  from  curves  presented  in 
Reference  4  for  low  speed  but  supercritical  Reynolds  numbers  as  a  basic  set  for  this 
investigation.  In  this  section  we  discus  an  alternative  approach  to  the  numerical  approach 
employed  in  the  Fortan  program  of  Appendix  A,  which  may  not  necessarily  be  as  precise  as 
that  of  the  code,  but  which  can  be  advantageous  for  engineering  applications.  The  data  as 
read  and  smoothed  data  are  given  above.  Some  of  the  points  as  read  are  plotted  in  figure 
B.l .  A  curve  has  been  fitted  to  these  data  in  order  to  give  them  an  approximate  analytical 
representation  over  the  limited  range  of  importance  for  inception  in  order  that  the  analytical 
approximation  can  be  a  reasonably  close  approximation.  The  region  of  constant  pressure 
indicates  roughly  the  position  of  the  laminar  separation  bubble  on  the  cylinder  which  is 
followed  by  boundary  layer  separation  and  a  turbulent  free  shear  layer  bounding  the  wake 
behind  the  cylinder.  Figure  B.  1  appears  to  span  the  range  of  arc  lengths  of  interest  for 
possible  vaporous  cavitation  bubble  growth. 

The  derivation  of  a  fencing  function  requires  that  the  pressure  distribution  be 
convened  from  the  dimensionless  arc  length  on  the  body  to  the  dimensionless  laboratory 
time  t,,  as  seen  by  a  typical  bubble  moving  with  the  vortex  sheet  which  approximates  the 
boundary  layer  in  this  study.  We  can  write  the  laboratory  time  as  in  Equation  (2.15): 

s 

■*=  J  7r%r  •  (BI) 

1.15 

were  the  lower  limit  denotes  the  first  abscissa  of  the  seventeen  points  in  Figure  B.l.  The 
upper  limit  refers  to  each  successive  point  Rather  than  use  the  curve  fit  for  the  pressure 
distribution,  we  have  used  the  data  as  shown  by  the  points  in  Figure  B 1  and  the  trapezoidal 
rule  in  order  to  correlate  all  plotted  points  with  tg.  The  seventeen  points  are  shown  in  the 
plot  of  figure  B2  in  which  tg  is  the  abscissa.  This  correlation  was  fitted  by  a  polynomial  as 
also  shown  in  die  figure.  From  these  correlated  (s,tg)  points,  each  Cp(s)  point  was  assigned 


Pressure  coefficient,  Cp(s) 
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Figure  B 1  Pressure  Distribution  Data  and  Approximating  Curve  Fit  for 
a  Circular  Cylinder  at  Supercritical  Reynolds  Number. 


Arc  length  on  cylinder,  s  radians. 
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Figure  B2  Correlation  of  Arc  Length  to  Dimensionless 
Laboratory  71036  for  a  Circular  Cylinder. 
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a  value  of  tg  in  order  to  construct  a  CpCtj)  plot  as  shown  in  figure  B3  .  Again,  an  approxi¬ 
mating  curve  fit  was  found  for  these  points  in  order  to  permit  an  analytical  representation  for 
the  forcing  function, 

F0s)=  }W,[-  Cp(ts)-K). 

As  explained  in  Section  2.2,  when  one  traces  the  pressure  distribution  starting  at 
the  stagnation  point,  the  forcing  function  has  its  origin  at  the  first  zero  of  the  quantity 
[Cpfts)  -  K].  The  laboratory  time  for  this  point  will  be  called  tg  and  it  can  be  found  from 
a  trial  and  error  solution  of  the  equation, 

Fi(ts)  =  »!  +  a2to  +  a3sin(a4to  +  a5)  -  K  =  0,  (B.2) 

in  which  the  notation,  F^tj),  indicates  a  radial  Weber  number  of  ~ ,  or  more  simply, 

the  quantity,  -  Cp  -  K.  The  solution  of  (B.2)  was  found  for  the  five  successive  K 
values  of  2.1, 2.2,  2.3, 2.4  and  2.5,  as  illustrated  in  Figure  B4.  A  cubic  fitting  curve 
for  tg(K)  is  also  shown  in  the  figure.  Then  we  can  put  tj  =  tg(K)  +t  in  Equation  (B.2) 
in  order  to  write  the  forcing  function  for  various  t  as 

Fj(f)  =  aj  +  a2(to  +  f)  +  a3sin(a4f +  a4tg+  a5)  -  K .  (B.3) 

A  short  Fortran  program  was  written  to  find  the  values  of  tg  and  evaluate  Fj(t)  for 
the  five  K  values.  Figure  B5  shows  the  function  F;(f)  although  the  abscissa  is  labeled  as 
tg  in  order  to  emphasize  the  fact  that  the  forcing  function  is  measured  in  units  of  laboratory 
time.  Figure  B.5  also  shows  the  start  of  laminar  separation  at  the  junction  of  the  curves 
from  equation  (B.3)  and  the  straight  segments,  having  ordinates  at  2.225  -  K.  At  any  K  this 
ordinate  is  constant  because  once  a  nucleus  finds  itself  in  the  laminar  separation  bubble  it 
will  grow  in  place  by  air  diffusion  until  it  becomes  large  enough  to  be  swept  into  the 
turbulent  free  shear  layer  downstream  of  the  cylinder.  This  approximate  laminar  separation 
point  is  defined  in  terms  of  the  shifted  laboratory  time  by 

?s  =  0.365  -to(K). 


(B.4) 


Pressure  coefficient,  Cn(t~) 
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0.0  0.1  0.2  0.3  0.4  0.5 


Dimensionless  laboratory  time,  t  s . 


Figure  B3  Pressure  Distribution  on  the  Cylinder  in  Terms  of 
Dimensionless  Laboratory  Time. 


Laboratory  time  ,  t0,  at  which  F(ts)  = 


36 


Figure  B4  Abscissa  of  Forcing  Function  Origin  in  Laboratory  Time 
for  a  Range  of  Cavitation  Numbers. 
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Figure  B5.  The  Bubble  Forcing  Function  for  Various  Cavitation  Numbers. 
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A  further  calculation  is  required  in  order  to  evaluate  the  initial  condition  u^O)  =  C  in 
shifted  time  coordinates  in  accordance  with  equation  (2.46).  The  evaluation  depends  upon 
the  derivative  dCp/ds  at  the  point  Sq  =  SofoCK)].  The  calculation  is  facilitated  by  use  of 
the  additional  fitting  function, 

s0(t0)  =  1.478  +  1.8797to  -  0.073  ltj.  (B.5) 


For  then  at  any  value  of  K  the  derivative  dCp/ds  can  be  evaluated  from  the  curve  fit  of 
Figure  B.l  at  the  correct  value  of  % 

In  the  special  case  when  uv  =  1,  we  have  1+  q  =  y  and  Equation  (2.46)  can  be 
written  as 


u1(0)  =  C  =  ^WrVY(l+K) 


dCp 

L"  ds . 


(B.6) 


dui 


Cp  =  -K 


Equation  (B.6)  pertains  to  ~r~  evaluated  at  tg  s=  to  (or  1  =  0  in  shifted  coordinates.) 

Qts 

When  uv  =  1  these  results  can  be  used  to  obtain  a  closed-form  approximate  solution  for 
bubble  growth  in  this  special  case. 

Or  more  generally,  equations  (B.3),  (B.4)  and  (B.5)  can  be  used  to  replace  some 

of  the  purely  numerical  parts  of  a  code  such  as  that  discussed  in  Appendix  A.  In  such  a 
procedure,  we  recall  that  since  Fj  holds  for  a  Weber  number  such  that  ( jWr )  =  1  the 


forcing  function  is  calculated  from 


R>s)  =  J  WrF,(ts).  (B.7) 

Morover,  although  the  present  example  supposes  that  a  laminar  separation  bubble  is 
present,  there  appears  to  be  no  reason  to  restrict  the  preceding  approach  to  shuch  flows.  The 
method  should  apply  to  unseparated  boundary  layer  flows  also. 
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APPENDIX  C 

LIMITING  VALUE  OF  THE  OSCILLATION-PERIOD  PARAMETER  X(uv) 


A  convenient  starting  point  for  the  present  analysis  is  equation  (3.20)  which 
defines  the  oscillation  period  parameter  as 

M»»>  -  =  2Kuv),  (3.20) 

where  the  half-period  integral,  I(uv),  is  defined  by  equation  (3.15)  as 

(3.15) 

l 

We  seek  the  limit  of  I(uv)  as  uv  — >  1.  Consequently,  as  we  have  seen,  um  =  1  +  ym  and 
£  =  1  +  y,  where  um  is  the  amplitude  of  oscillation  and  0  <  y  <  ym  <  1  because  the 
amplitudes  are  small  when  uv  is  in  the  neighborhood  of  unity.  Therefore  we  can  rewrite 
equation  (3.15)  as 
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b0=  -  +  '\](-)-V-.  Therefore,  in  ascending  order,  the  three  roots  are  0,<  ym,<  b. 
Next  we  can  put  y  =  ym  4  in  order  to  write  (C.l)  as 


_ (i  +  Ym  £)3 _ ^ 

(2/3)uXW-5>[i-(y«,VM 


(C.2) 


Now  we  expand  the  factor, 

— —  +  y™^  -  =l+~(3  +  “)y  £  +  -  (1  +  —  +  —  )y2  ?2  + 

ri  K/u  \H/5  1  +  9  V-5  f  u  jymS  +  8'1  +  K.  +  .  2  +  " 


[l-Cy^/bo)]'/2 


8V  bo  b2"m 


which  enables  us  to  decompose  (C.2)  into  three  integrals.  These  become 


(2/3)  u>0 


11  +  2  (3  +  b0  )ymI2  +  \  0  +  b0  +  b2  )ymJ3 


(C.3) 


in  which 


l 

/* 


1 

r 


and  I3  = 


V50-5)' 


Clearly,  Ij  =  7t  and  I2  =  -  .  The  third  integral,  I3,  can  be  dealt  with  if  the  radicand  in 

l  1 

its  integrand  is  written  as  ~  -  (q  -  ~)z.  Then  the  numerator  becomes 

Therefore  I3  can  be  expressed  as  the  sum  of  three  integrals,  which  upon  integration  give 

_  L/2E  ZE\  5.  _  3£ 

l3-8(2  +  2'  +  4“8  ‘ 

These  values  of  Ij,  I2  and  I3  can  now  be  put  into  equation  (C.3).  The  result  is 


(2/3)  u^b0 


i  +  1i(3+^)ym+ 


(C.4) 
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As  uv  ->  1,  k  -*  3,  |i.  ->  0,  ym->  0  and  b0-»  3.  Therefore  I  -»  .  From  this  result  and 

equation  (3.20)  it  follows  that 

X(l)  =  jcV2  =  2.2214415  .  (C.5) 

Next  we  repeat  the  analysis  using  a  different  method.  Therefore  we  return  to  the 
integral  I  and  write  it  as 


ym 


0  +  y) 


y(ym-y)(bo-y) 


dy 


-iVIW  <c6) 


where 


ym 


II- 


V: 


■Q  +  y) 


y(ym-y)(b0-y) 


— r  dy  and  I2  = 


y**1  ___ _____ 

\j  —  y(i  +  y)~ 
J  V(ym-y)(b0- 


y) 


dy 


0  0 

are  complete  elliptic  integrals.  The  properties  of  such  integrals  are  discussed  in  detail  by 
Byrd  and  Friedman1.  In  terms  of  the  roots  displayed  by  the  factors  in  the  integrands  we 
can  write 


k2  =  Mko±H 

b0(ym+l)’g 


,  == . — ,  a2  =  “f-  and  <j> 

Vb0(ym  +  1)  bo 


•  /  b0(ym-y) 

Sin  V  ym(bo-y) 


For  complete  integrals  the  lower  limit  of  integration  is  y  =  0  and  <j>  but  if  y  =  ym 

*  =  0.ThenwchaveI1  =  2-^^^fF(|,  k)  +  (%,“ '+*})  II(f  .a^.k) 

where  F(^,  k)  is  the  complete  elliptic  integral  of  the  first  kind  and  where  IT(  j  ,a2,k) 

is  the  complete  elliptic  integral  of  the  third  kind.  Then  we  take  the  limit  as  uv— >  1, 
ym->  0,  k  -*  0,  and  a  ->  0.  Hence  F  (|,  0)  =  |  and  II(§,  0,0)  =  F(j,0)=j  . 

Recalling  that  b0=  3  in  the  limit,  we  see  that 


1Paul  F.  Byrd  and  Morris  D.  Friedman,  Handbook  of  Elliptic  Integrals  for  Engineers  and  Physicists, 
Springer  Verlag,  Berlin,  19S4.  See  p.  1 16,  *255.20  and  *255.1 1. 
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r  2  It  _v  2lt 

I,‘  V52  (I  ~  3)_“  V3  • 

From  Byrd  and  Friedman  we  find  that  I2  can  be  expressed  in  terms  of  elliptic 
integrals  of  the  first,  second  and  third  kinds.  After  some  manipulation  they  can  be 
written  in  the  form  of  complete  integrals  as 

h  =  2  V(ym+l)bo  ' j-Etfjc)  +  (1  +  ^~-)F(f*)  +  (1  -  jj^y)II(f,a2,  k)  j. 
In  view  of  the  limiting  values  of  the  several  parameters  one  finds  that 


I,  =  2VT 


(  n  5k  3k\ 

[  2+  2  “  2  J" 


tr/3 


Then  from  equation  (C.6)  with  uv  =  1  we  find  that 


1  = 


+  K  V3 


K 


in  agreement  with  the  value  of  I  given  above.  Therefore  the  value  of  I  given  in  equation 
(C.5)  is  confirmed. 
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APPENDIX  D 
BACKGROUND  NOTES 


D  1  Basic  Bubble  Dynamics 

The  well  known  Rayleigh-Plesset  differential  equation  governs  the  growth  or 
collapse  of  a  spherical  vapor  bubble  in  a  perfect  fluid  when  the  static  pressure  surrounding 
the  bubble  varies  with  the  time,  t,  [16].  In  the  present  investigation  we  shall  consider  the 
vaporous  growth  of  such  bubbles  in  a  flow  in  which  the  the  bubble,  in  moving  along  the 
surface  of  a  submerged  body,  is  influenced  by  its  pressure  distribution.  The  bubble  will 
contain  air  and  water  vapor  so  that  the  dynamic  force  balance  is 

p(RR  +  §R2)=Pi-2S-p(t).  (d  i) 

In  this  equation  R(t)  is  the  bubble  radius,  p  is  the  liquid  density,  a  is  the  surface  tension 
and  p(t)  ir  the  static  pressure  distribution  on  the  body  as  seen  by  the  moving  bubble.  The 
pressure  inside  the  bubble  is  pj  and  if  the  bubble  is '  jothermal ,  Pi  =  pv  +  Pa  (Ro/R)3.  Here 
pv  is  the  vapor  pressure  and  pa  is  the  air  pressure  in  the  bubble  at  radius  Rq  in  the  free 
stream.  The  radius  Rq  is  the  size  of  a  typical  microscopic  cavitation  nucleus  and  pa  is 
related  to  the  dissolved  air  content  in  the  water.  Therefore  the  basic  equation  of  present 
interest  is 


p(RR  +2  R  )  =  Pv  +  PaOV^)  -  2o/R  -  p(t).  (D2) 

Next  we  shall  introduce  dimensionless  variables.  We  start  by  letting  R  =  r  Rq  and  observe 
that 

P(0  =  Po  •  I  Pvo  Cp  “d  *^at  Pv  =  Po  +  2  Pvo  K  *  where  Cp  is  the  pressure  coefficient  on 

the  submerged  body  and  K  is  the  cavitation  number.  Therefore  we  have  py  -  p(t)  =  - 
2  PV0  (Cp+  K).  Substitution  of  these  definitions  into  the  differential  equation  leads  to 
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pRo2 


2a/R0  fr?+lf2)-W?'7  +  iw^cr^ 

pv? 


(D3) 


where  the  definition  of  Wr  is  Wr  =  —7= —  ,  the  radial  Weber  number.  We  shall  define  the 

O/Kq 

Pa  PR°2 

air  content  parameter  by  y=  2.0 /Rq  ^  also  observe  that  the  factor  2aiR^  ^as  t^ie 

dimension  of  [Time]2.  Therefore  we  shall  define  a  new  dimensionless  time  x  by  x  = 
t  <\j  (2a)/(p  R^) .  Consequently  we  can  write  the  dimensionless  equation  of  motion  as 


tM-  +3(*)2,x  .  I  +F(t) 

r  dl*  2  W  r3  r  '  >• 
where  the  forcing  function  is  defined  by 

F(t)  =  4  Wr  [-  Cp(t)  -  K]. 


(D  4) 


(D5) 


D2  Laminar  Separation  and  Cavitation  on  Circular  CvUnoers 
Cavitation  inception  observations  on  circular  cylinders  have  been  reported  by  Ihara  and 
Murai  [3].  Their  experiments  were  designed  to  see  whether  or  not  laminar  separation  plays 
a  role  in  the  onset  of  cavitation  cm  smooth  cylinders  in  critical  and  supercritical  flow  ranges 
and  on  cylinders  having  a  boundary  layer  trip.  They  observed  inception  resembling  bubble¬ 
ring  cavitation  in  the  reattachment  region  when  a  laminar  bubble  was  present  at  Reynolds  as 
high  as  3.27x10s  cm  a  smooth  cylinder.  They  called  this  cavitation  “bubble-line”cavitation. 
Ordinarily,  supercridal  flow  is  encountered  at  Reynolds  numbers  of  3x10s  or  greater.  No 
tunnel  wall  effect  corrections  are  noted  in  the  paper.  Such  effects  must  have  been  present 
because  the  tunnel  height  was  5.42  cylinder  diameters,  and  at  a  Reynolds  number  of 
3.3x10s  a  minimum  pressure  coefficient  of  -2.9  is  reported.  This  Cp  magnitude  is 
somewhat  larger  that  those  ordinarily  observed  in  the  absence  of  wall  effects.  This  fact  is  of 
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minor  importance  because  these  authors'  measurements  do  establish  the  role  of  laminar 
separation  on  smooth  circular  cylinders  in  this  Reynolds  number  range;  but  it  suggests  to 
us  that  we  should  use  the  pressure  measurements  of  Gowen  and  Perkins  [3]  which  appear 
to  be  free  of  wall  effects. 


